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Abstract 



A new method for controlling harmonic generation, in the framework of quantum optimal 
control theory (QOCT), is developed. The problem is formulated in the frequency domain 
using a new maximization functional. The relaxation method is used as the optimization 
procedure. The new formulation is generalized to other control problems with requirements 
in the frequency domain. The method is applied to several simple problems. The results 
are analysed and discussed. General conclusions on harmonic generation mechanisms are 
obtained. 
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Preface 

General remark: the eigenstates of the Hamiltonian, in the various problems, will be 
denoted by |(^n); the index n represents the eigenstate with the eigenenergy En, where: 

Ei < Ei+i z = 0, 1, . . . 

General remarks for all numerical results: 

• Atomic units are used throughout. 

• The important details of the problem and the computational process are presented in 
a table. Most of the notations are defined in the text. The meaning of all notations 
is described in the following table: 
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Notation 


Description 


Ho 


the unperturbed Hamiltonian 




the dipole moment operator 


iV'o) 


the initial state vector 




the initial wave function in the x domain 


10) 


the target state vector 


(j){x) 


the target wave function in the x domain 


T 


the final time 


a 


the penalty factor 


fM 


the filter function of the forcing field 




the filter function of the dipole moment expectation value 


e\t) 


the initial guess of the field, in the time domain 




the initial guess of the field, in the frequency domain 


L 


n = L is the index of the maximal allowed eigenstate 


In 


the penalty factor of the forbidden state 


K 


the penalty factor of [ '^'''^]f^ ) 




the initial guess of K, for the relaxation method 


X domain 


the domain of the x grid 




the number of points in the x grid 


tolerance 


the tolerance of the convergence of the field (see App. B) 



u{x) denotes the Heaviside step function: 




X < 
< X 
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Chapter 1 
Introduction 



Harmonic generation is a process in which the frequency of an electromagnetic radiation 
is multiphed by a quantum system. This phenomenon occurs when the system emits 
radiation at a frequency higher than that of the incident radiation. 

In the harmonic generation process, the incident forcing electromagnetic field produces 
oscillations in the dipole moment expectation value {fi){t) of the system. According to 
Maxwell equations, an accelerated charge emits electromagnetic radiation, proportional 
to the acceleration of the charge. Hence, an oscillating quantum system emits radiation, 
proportional to the acceleration of {fi){t). The spectrum of the field emitted by the system 
consists of the same frequencies as the spectrum of the {fi){t) oscillations. Harmonic 
generation occurs when the {fi){t) spectrum contains frequencies higher than that of the 
forcing field. 

Typically, the most important component of the (fi) (t) spectrum is the frequency of the 
forcing field. This component represents the linear response of the system to the radiation. 
The harmonic generation phenomenon originates from non-linear effects. Typically, these 
are small in magnitude compared to the linear response. 

The harmonic generation phenomenon has been utilized since the invention of the laser 
to convert radiation from a laser source to a frequency higher than that available from this 
source. Most commonly, the laser radiation frequency is multiplied by a factor of 2 (second 
harmonic generation) or 3 (third harmonic generation). 

In the late 1980's, a new harmonic generation phenomenon was discovered. In this 
phenomenon the incident laser radiation is multiplied by factors of tens, even hundreds. 
This phenomenon is known as "high-harmonic generation". The high-harmonic generation 
phenomenon is the key element in the production of attosecond laser pulses; these are of 
extreme importance for investigation and control of electronic processes [15]. 

As mentioned, the harmonic generation phenomenon is small in magnitude. In addition, 
the harmonic generation spectrum often consists of many frequencies, while interest is in 
only a single frequency, or a region in the spectrum. It is of great interest to control the 
harmonic generation phenomenon in order to increase the intensity of the emitted field at 
the frequencies of interest. 

The problem may be addressed by means of the quantum control discipline. It deals 
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with the task of controlhng quantum systems by an electromagnetic radiation using de- 
signed field sequence shapes. The control of the harmonic generation phenomenon may 
be achieved by designing an appropriate time-shape of the incident laser field. This field 
sequence will extend over a spectrum of relatively low frequencies instead of being limited 
to a monochromatic radiation. 

In order to control a quantum system it is desirable to find the optimal field for the 
problem of interest. There are two approaches for seeking the optimal field for quantum 
control problems: 

1. Experimentally, by a sophisticated trial and error process using a genetic algorithm 
(see [1]); 

2. By theoretical calculation, using our knowledge about the quantum system. 

In the experimental approach all the aspects of the problem are taken into account, unlike 
in the theoretical approach. This makes the experimental approach much more accurate 
and efficient. The advantage of the theoretical approach is the possibility of investigating 
the control mechanism. This is the main importance for providing a theoretical method 
of calculation. The theoretical calculation may also provide a good starting point for the 
experimental genetic algorithm search. Frequently, a good starting point is necessary for 
the success of this method. 

Quantum optimal control theory (OCT) is the most successful theoretical method avail- 
able for finding an optimal field in quantum control problems. It is sometimes shortened as 
QOCT. In the framework of QOCT, the problem is formulated as a maximization problem 
using the calculus of variation formalism. 

Great progress in the task of controlling high harmonic generation has been achieved 
in the last decade by using the experimental approach [15]. However, a theoretical method 
of calculation for harmonic generation is still missing. The main purpose of the present 
work is to fill this gap. 

The goals of our research are: 

1. Developing a procedure for finding an optimal field for harmonic generation, in the 
framework of QOCT; 

2. Finding a way for dealing with control problems with frequency requirements, in the 
framework of QOCT; 

3. Exploring new mechanisms of harmonic generation. 

In Ch. 2, we present the relevant theoretical background for the present work. In Ch. 3, 
the new method is developed. In Ch. 4, the new method is applied to several harmonic 
generation problems. The results are discussed, and general conclusions on mechanisms of 
harmonic generation are obtained. 
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Chapter 2 
Background 



In this chapter we present the basic background in QOCT, needed to understand the 
present work. We also present existing works that are related to harmonic generation task. 
We can divide our task into two distinct parts: 

1. Imposing a restriction on the spectrum of the forcing field 

2. Controlling the spectrum of the oscillating dipole, which is the spectrum of the emit- 



In this chapter, and also in the following one, we shall treat these two parts separately. 

In sections 2.1-2.3, we present the basic relevant background in QOCT. Section 2.4 deals 
with the existing works related to the first part of our task, mentioned above. Section 2.5 
deals with a work related to the second part, and the failure of an attempt to combine the 
two parts for harmonic generation control. 

2.1 OCT of a target operator in the final time 

The simplest and most common task in QOCT, is seeking a time dependent forcing field 
that maximizes the expectation value of an operator in a final time (see [1, 2, 6, 7], in 
more detail). Hence, the formulation of this problem is the natural starting point for any 
discussion in QOCT. 

Consider a quantum system in an initial state: |^(0)) = {ipo), with the unperturbed 
Hamiltonian: Hq, under a forcing field: e(t); the OCT formulation translates the above 
mentioned maximization requirement into the maximization of the following functional: 



ted field 




(2.1) 



where O is an arbitrary operator, and T is the final time. \ip{T)) depends on e{t) in a rather 
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complicated way, through the Schrodinger equation, under the given initial condition: 



= -^H(t) m)) , 1^(0)) = \^o) (2.2) 

H{t) = Ho - Kt) 

where fi is the dipole moment operator. (Atomic units are used throughout, so we set: 
h = 1.) We have used the dipole approximation for H(t). 

The most common target operator is the projection operator on a given target state, 
which will be denoted as |0): 

P0 = 10) (01 (2.3) 

In this case, Jmax will attain its maximal value when: \ip{T)) = e*^ |0) [6 is an arbitrary 
phase). This kind of target operator is used when we want to find an optimal field for a 
state-to-state transition. 

If we try to seek a maximum of Jmax, "we obtain non-physical fields of very large, short 
timed pulses; the reason is that the solution of this maximization problem is singular, 
with an infinite field. To get a physical solution, we must restrict the intensity of the 
field. It is also desirable, practically, to use lasers of as small an intensity as possible. 
This requirement can be achieved, by adding a "penalty" term to the functional object of 
maximization: ^ 

JpenaiW)] = -a e\t) dt a>0 (2.4) 

Jo 

This term will be maximized when e(t) is minimal. It "penalizes" the object of maximiza- 
tion for using high intensity fields, a is a "penalty factor" - a positive constant, whose 
value determines the "cost" of large intensities. The suitable value for the problem usually 
has to be determined by a trial and error process. 
Let us define the functional: 

JmpHt)] = JmaxHt)] + JpenalHt)] (2.5) 

Our problem is to find a function e{t) that maximizes Jmp- This kind of a problem belongs 
to the category of variational problems. It can be handled most conveniently by the 
functional derivative formalism. The condition for an extremal is: 

(2,6) 

The complicated explicit dependence of Jmax on e{t), makes dealing with Eq. (2.6) rather 
inconvenient (see [12] for such an approach). The most convenient way to handle the prob- 
lem is by using the Lagrange-multiplier method; it enables us to treat implicitly dependent 
variables as independent, in the first stage. Eq. (2.2) will be treated as a constraint, which 
will be enforced later. 

Before we proceed, we have to choose the variables that will be treated as independent. 
We should notice, that |'0(t)) is complex, and is composed of two distinct variables: the real 
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and imaginary parts. Moreover, the constraint equation, (2.2), consists of two conjugate 
constraint equations — the equivalence relations between the real and imaginary parts of 
the expressions on both sides of the equation. There are also two initial conditions. We 
could have treated 'Re\ilj{t)) and Im|-?/;(t)) as our independent variables; nevertheless, it 
appears to be much more convenient to treat |'0(t)) and {ip{t)\ as independent. Now, we 
can use Eq. (2.2) as is, as our first constraint equation and initial condition; the additional 
constraint equation and initial condition is given by the complex conjugate of Eq. (2.2): 



d{m\ 

dt 



(2.7) 



(2.7) ensures that (^^(t)! will always be the complex conjugate of lipit)), as required. 

Now, we modify Jmp by adding a constraint functional term. First, we write equa- 
tions (2.2), (2.7) in the following way: 



d\m) 

dt 

d{m\ 

dt 



-lUit) \m) = o 



(2.8) 
(2.9) 



Note that equations (2.8), (2.9) refer to all t in the interval: < t < T; hence, they 
impose distinct constraints on the state at all time points of the interval. According to 
the Lagrange-multiplier method, we have to add to the functional a distinct term for every 
constraint. Let us start with the constraints imposed by (2.8): We add to Jmp a continuous 
summation over all the LHS expressions of the constraints in all t, each multiplied by its 
own Lagrange-multiplier, (x(^)l; the sequence of Lagrange-multipliers forms a continuous 
function of t in the interval (see [1, Ch. 2]). The additional term for the constraint in (2.8) 
takes the form: _ 

i+^H(t) 



T 



x{t) 



dt 



i){t) ) dt 



(2.10) 



We take (x(^)l to be a complex function, like a wave function. |x(^)) is sometimes called: 
"the conjugate function" of \ip{t)). 

We treat the additional term for the constraint in (2.9) in the same way. We have 
to use another Lagrange-multiplier function of t for this term; since we took (x(^)| to be 
complex, we may treat |x(^)) as an independent function. The additional term is: 



^ ' 'I - zH(t) ) m 



xit) ) dt 



(2.11) 



Note that (2.11) is the complex conjugate of (2.10). 

The overall additional constraint term for the functional object is: 



J con 2Re 



T 



x(t) 



i){t) ) dt 



(2.12) 
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The overall new object of maximization is the following functional: 



J — Jmax ~l~ Ji 



penal 



Jr. 



i){T) O il){T))-a / e^{t)dt-2Re 



T 



T 



x{t) 



d_ 
di 



zH(t) 



'0 ^0 

The conditions for an extremal are given by the following equations: 

6J 







6eit) 
6J 

mrY) 

6J 













^{t) ) dt (2.13) 

(2.14) 
(2.15) 
(2.16) 
(2.17) 
(2.18) 



6{i;{T)\ 

together with the constraints, (2.2) and (2.7). All functional derivatives are taken while 
holding all other variables constant. 

This set of equations forms the basis for the derivation of the so called: "Euler-Lagrange 
equations" of the problem . The resulting equations are (see [1, 2] for example, for more 
details) : 



dt 
d\x{t)) 

dt 
H(t) = 

e(t) = - 



= -zH(t) m)) , 

= -tU{t) \x{t)) , 

Ho - fie{t) 
lm{x{t)\fi\iit)) 



\m) = m 

|X(T)) = 6|^(T)) 



a 



(2.19) 
(2.20) 

(2.21) 



A solution e{t) which satisfies all these equations is an extremal field, which gives a maxi- 
mum for the functional J, and hence, also for Jmp (a minimum is almost never encountered 
in this kind of problems). 

There is no general method to solve this set of equations analytically. Hence, we have 
to employ numerical methods to find a solution. Several methods will be presented in 
section 2.3. 



2.2 OCT of a time dependent target operator 

The task that was described in the last section refers to the control of the system state in 
a single time point: t = T. A more difficult problem is to control the state of the system 
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over an interval of time. This is the kind of problem that we deal with in the present work. 
In this section, we describe the existing method for the control of such kind of problems 
(see [14, 2] for more details). 

The object of maximization is the expectation value of a time dependent operator, 
during the time interval: < t < T. The Jmax from Eq. (2.1) is replaced by (see [2]): 



J rr 



T 



ip{t) 6{t) ip{t)) w{t) dt 



(2.22) 



where w{t) is a weight function, which satisfies: 



w{t) dt 



w(t) determines the relative importance of the maximization target during different parts 
of the time interval. The functional will be maximized when the integrand is maximized 
at all time points in the interval. 

The method is intended for a positive-semi definite 0(t), in which a greater expectation 
value indicates greater success. This is not always the case — it may occur that we want 
the expectation value to vary in a predefined path, and not to be maximized at all t. In 
such this method is not helpful. 

A common example for a time dependent target operator is a time dependent projection 
operator: 

p^it) = \m{m (2.23) 

This target operator should be used when we want to force the system to follow a predefined 
path, defined by the sequence of wave-functions: |0(t)). 

The other parts of the object of maximization J, i. e. Jpenai and Jcon, are the same as 
in Eqs. (2.4), (2.12). The overall object of maximization is the following functional: 



J — Jmax ~l~ Jpenal ~l~ Jcon 

{^^(t) 6{t) i){t)^w{t)dt-a e^{t)dt-2ReJ^ ( x(t) 



d_ 

dt 



iU{t) 



i){t)j dt 
(2.24) 

The conditions for an extremal are the same as in Eqs. (2.14)-(2.16), together with 
the constraints: (2.2) and (2.7). The resulting Euler-Lagrange equations are (see [2] for 
details): 

dm)) 



dt 
d\x{t)) 
dt 



-zH(t) m)) , 
-tYi{t)\x{t))-w{t)6{t)\m) 



\m) = m 

\X{T))=0 



Hit) = Ho - fie{t) 

f.^ _ Im(x(t) \fi\ip{t)) 
'^^^ ~ a 



(2.25) 
(2.26) 

(2.27) 
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Note that (2.25) and (2.27) are the same as (2.19) and (2.21). However, (2.26) is different 
from (2.20); we have here an additional, inhomogeneous, \ip{t)) dependent term on the 
RHS of the equation. This form of equation is known as "the inhomogeneous Schrodinger 
equation". The boundary condition is also different. It does not come from the conditions 
for maximization in Eq. (2.14)-(2.16); this boundary condition represents the "natural 
boundary conditions" of the problem (see, for example, [4, Ch. 2, Sec. 3]). The fact 
that the Lagrange- multiplier |x(^)) vanishes originates from the independence of Jmax 
on \tlj{T)) (more precisely, the dependence is infinitesimal, and does not contribute to the 
boundary condition in T). This should be compared with the boundary condition in (2.20), 
where Jmax does depend on \ip(T)). 

This boundary condition is problematic because it imposes on the field the condition 
of vanishing at the final time: 

6(T) = 

Our experience shows that this condition is the source of difficulties in achieving control 
in the neighbourhood of t = T. 

Another version of the formulation of the problem is available in the literature ([4, 
Ch. 2], [S]), in which this difficulty is eliminated; J is modified by adding the following 
boundary term: 

Jbound = U{T) 6{T) i,{T)) > (2.28) 



K is a positive parameter, that determines the relative importance of the boundary term. 
The modified object of optimization is: 

J — Jmax ~\~ Jbound ~l~ Jpenal ~l~ Jcon (2.29) 

Now we have additional conditions for extremum: (2.17) and (2.18). 

The resulting Euler-Lagrange equations are the same as (2.25)-(2.27), except the bound- 
ary condition in (2.26), which is replaced by: 

\X{T)) = «:6(T) \^{T)) (2.30) 

The disadvantage of this approach is that it gives exaggerated importance to the target 
at T, compared to the other time points. 

2.3 Numerical methods for the maximization of the 
functional 

In this section, we present a number of numerical methods available for the maximization 
of the J functionals mentioned in the two previous sections. 
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2.3.1 The naive approach 

As a first thought, we may propose the iterative scheme, described in the following steps 
[1, 7]: 

1. Guess a field sequence: e(t). 

2. Repeat the following steps, until convergence: 

(a) Propagate \ip(t)) forward from t = to t = T, using (2.19), with e(t). 

(b) Set according to the boundary condition in (2.20), (2.26) or (2.30). 

(c) Propagate |x('^)) backward from t = T to t = 0, using (2.20) or (2.26), with e{t) 
(and Itpit)) from step 2a, for (2.26)). 

(d) Update e(t) according to (2.21), using the sequences of I'lpit)) and |x(^)) from 
steps 2a, 2c. 

Unfortunately, this simple scheme seldom converges. More elaborate methods are 
needed for control problems. 

2.3.2 The gradient methods 

The gradient methods are a family of general methods of optimization; they use information 
about the derivatives of the object of minimization/maximization with respect to the 
variable to be optimized (see [4, 5]). In our case, we talk about functional derivatives 
of J with respect to the control field e(t). This information is used to "climb up" in the 
hypersurface of J vs. e{t) (at all time points). 

The simplest gradient method is the so called: "first-order gradient method" , or "steep- 
est descent / accent method" . As its name indicates it consists of the information from the 
first order derivative, or the gradient, of the object of maximization. In our case, the 
gradient is: 

V.w J =j-^ = -2 [ae{t) + lm{x{t) ^(t))] (2.31) 

The gradient is a vector (in our case, a continuous vector, i. e. a function) that points in 
the direction of the maximal increase of the object of maximization. The method is based 
on following this direction. In our case, we update e(t) according to the following rule [6]: 

K >0 (2.32) 

where i^' is a positive parameter, which determines the rate of propagation in the direction 
specified by Ve(t)J- We repeat the process until convergence. 

K must be small enough for the first order approximation to be satisfactory; if K is too 
large, the sign of the gradient might change in the way to the new e{t), and the value of J 
will not necessarily increase. On the other hand, if K is too small, the rate of convergence 



e-<="'(t) = e°''^(t) + J^Ve(t)J 
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will be slow. An optimal value of K has to be found by a trial and error process. In 
addition, it should be varied during the process of optimization, from greater values in 
the early iterations, to smaller values close to the maximum. A commonly used method 
for searching for an optimal value of K is called "line search" . The line search has to be 
reimplemented in every iteration. 

Let k denote the iteration index; the whole procedure is summarized in the following 
scheme: 

1. Guess a field sequence: e'^°)(t). 

2. Propagate [^/''•^^(t)) forward from t = to t = T, with e'^^\t). 

3. Calculate J(°) with \i)^'^\t)) and e(°)(t) (note that Jcon = 0, since the constraint is 
satisfied) . 

4. (k = 0) 

5. Repeat the following steps, until convergence: 

(a) Set |x''^H^)) according to the suitable boundary condition, using |'?/;(^^(T)). 

(b) Propagate \x^^H'^)} backward from t = T to t = 0, with e^''''{t). 

(c) Perform a line search, to find an optimal value for K; the line search involves J 
evaluations for various values of K, to be compared with J^^-*; these require the 
following steps: 

i. Set a new field, using Eq. (2.31): 



(2.33) 



ii. Propagate | (t)) forward from t = to t = T, with e*""'(t). 

iii. Calculate J*""^' with | (t)^ and e*""'(t), and compare with J^'^^ 

(d) When an optimal K was found, update all the variables according to this K: 

= e*"»'(t) |^('=+^)(t)> = |/""'(t)) Z'^+i) = J*""' (2.34) 

(e) (k = k + 1) 

The first order gradient method usually shows considerable improvement in the first 
iterations. However, it is known to have a very slow convergence rate when getting close 
to the maximum. Hence, we should avoid using it if possible. 

The "second order gradient method" , or "Newton method" , uses second order derivative 
information about the object of maximization, in addition to the gradient. It has good 
convergence characteristics near the maximum. This method is not suitable for QOCT 
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problems, for a reason that will be explained later; however, approximate versions of this 
method — "quasi- Newton methods" — can be employed successfully [11]. 

We start with the description of the regular Newton method for a nonlinear function 
of a single variable — f{x)- Suppose we have an initial guess xq, which is known to be 
close enough to the maximum; in "close" we mean that the negative sign of f"{x) does not 
change from the maximum to Xq. Now, we adjust a parabola to this point, which is the 
second order approximation for f{x): 



f{x) ^ p{x) = f{xo) + f'{xo){x - xo) + -f"{xo){x - Xof 
Then, we move to the x value of the parabola's maximum; it is easily found to be: 

x* = Xo- [fixoT^fixo) 



(2.35) 



We update the value of x to x*. We repeat the process, until convergence. 

If we deal with a function of a vector, the procedure is very similar; let F{\) be a 
function of the vector: 



Suppose we have an initial guess, vq, which is close enough to the maximum; in this case, 
the word "close" refers to the negative-definiteness of the Hessian matrix 5, which is defined 
by: 

We adjust a multidimensional paraboloid to the point vq, and we find its maximum. The 
expression for v at the maximum is very similar to (2.35): 



V* = vo - S-^VF 



(2.37) 



v=vo 



V is updated to v*, and the process is repeated until convergence. 

In our case, the Hessian is a "continuous matrix" , defined by a two variable function: 



s{t,t') 



5\J 



6e{t)6e{t') 



(2.38) 



The operation of the Hessian matrix on a vector in the discrete case, is replaced by the 
operation of an operator on a function in the continuous case; the Hessian operator is 
defined by: 



s{t,t')g{t') dt' 



(2.39) 
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where g{t) is an arbitrary function of t. The update rule is: 

(2.40) 

The main problem with the Newton-method is that often the Hessian is too complex 
to be computed easily. The expression for the first derivative in Eq. (2.31) is dependent on 
e{t) explicitly, and implicitly, through \x{t)) and \ijj{t)). The implicit dependence is very 
complicated; it follows, that the expression of the second derivatives, s{t,t'), is also very 
complicated. Hence, practically, the Hessian cannot be computed directly for all t, t', in 
every iteration.^ 

The quasi-Newton methods make use of an approximated Hessian. The most commonly 
used method is the BFGS method. It approximates the Hessian using an information from 
the gradient (see [5, 11] for more details). We will denote the approximated Hessian as: 

'-'ap- 

When using an approximate Hessian, Eq. (2.37) (in the discrete case) should not be 
used as is. Instead, we treat the vector —S'pVF as the direction of search for the new 
V. This vector has exactly the same role as the gradient vector in the first order gradient 
method. The update rule in the discrete case is: 

K>0 (2.41) 

v=vo 

and in the continuous case: 

K>0 (2.42) 

A line search is made to find an optimal value for K. 

The full process of the quasi-Newton methods, in the context of QOCT, is very similar 
to the scheme that was presented for the first order gradient method. The only difference 
is that we use another direction for search; Eq. (2.33) is replaced by the following: 

^trial ^ ^{k) _ g(fc)-l ^^^^ j| ik) ^2.43) 

where the Hessian approximation Sip\ and the gradient [Ve(f)J]^^\ are computed using: 

^You may ask: why do we treat |x(0) ^-nd as e(i) dependent when deahng with the second deriva- 

tives, while they were treated as t{t) independent when we dealt with the first derivative (in Eq. (2.31))? 
The answer is, that in fact, \x{t)) and \'^{t)) do depend on e{t); we ignore this dependence when handling 
with the first derivative, because this is exactly the role of the Lagrange-multipliers — they are adjusted 
in a way that makes Eqs. (2.15)-(2.18) true. Then, the full derivative of J with respect to e(i), coincides 
with the partial derivative in Eq. (2.31). However, when treating the second derivatives of J, this choice of 
the Lagrange-multipliers has no special significance, and the implicit dependence on e[t) cannot be ignored 
(see [.3]). 



e"--(t) = e°'^(t)-S-Ve(t)J 



16 



2.3.3 The Krotov method 



The most popular method for QOCT problems is called: "the Krotov method". It was 
first introduced, in the context of QOCT, in [9]. There are several variants of the method; 
here we present the variant that is used in the examples of this thesis [2]. 

The idea is not very different from the one of the naive approach; the difference between 
the methods is in the stage where the update of the field takes place. In the naive approach, 
the field is updated after the propagations of \ip{t)) and \xit)) were completed; in the 
Krotov method, the field is updated during the propagation, at every new time point, and 
the propagation proceeds using the new field. 

The procedure is presented in the following scheme: 

1. Guess a field sequence: e^°^(t). 

2. Propagate \4'^'^\t)'^ forward from t = to t = T, with e^°^(t). 

3. (k = 0) 

4. Repeat the following steps, until convergence: 

(a) Set according to the suitable boundary condition, using |^/;('^^(T)). 

(b) Propagate \x^''\t)'j backward from t = T to t = 0, with the field defined by: 

~(fc) ^ = \^ KJ\t-\y- \)/ 2.44 

a 

(c) Propagate \il}'^^^^\t)') forward from t = to t = T, with the new field, defined 
by: 

= Im(xW(t)|/x|^(^-+^)(t)) ^2^^^^ 
a 

(d) (k = k + 1) 

As can be seen, the procedure is considerably simpler than that of the gradient methods. 

The most important property of the Krotov algorithm is that it is monotonically con- 
vergent (with no need of a trial and error process, as in the gradient methods). The 
convergence analysis is presented in [10]. 

Contrary to the gradient methods, which are general methods, the Krotov method 
utilizes the special structure of the problem. The update rule of e{t) depends only on the 
functions in time t; this allows the update of e{t) during the propagations, in such a way 
that the resulting sequence of e{t) is consistent with the resulting sequence of the wave 
functions, and with itself. 

The rate of convergence in the Krotov method is known to be faster than that of 
the gradient methods, in the early iterations. However, when getting very close to the 
maximum, the convergence rate becomes slow. Hence, the Krotov method is problematic 
when a high fidelity solution is required. In [11] it was shown that the Krotov method 
degenerates into a first order gradient method close to the the maximum. 
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2.4 QOCT with a restriction on the field spectrum 



The task of restricting the spectrum of the forcing field has considerable importance in 
QOCT; the reason is that most of the computed fields turn out to be too oscillatory 
to be produced experimentally. Several methods have been proposed to overcome this 
problem [12]. Most of the methods are based on limiting the spectrum of the field. In this 
section we present the important methods. 

The methods were formulated in the context of the problem of Sec. 2.1. They can be 
used for other problems without additional mathematical complications. However, the rate 
of success of the method might not be the same. 

2.4.1 The method of Werschnik and Gross: Brute force 

Werschnik and Gross [2] use the Krotov method as a basis for the optimization process. 
However, the process is interrupted in every iteration by a spectral filtration of the resulting 
field. The filtered field is used in the next propagation. This amounts to replacing the 
update rule in (2.45) by the following: 



where T stands for the Fourier transform, for its inverse, lo is the frequency variable, 
and /(cu) is a filter function. 

This kind of interruption in the process is called "brute force" . The brute force method 
inserts into the process of optimization something that is not self-consistent with its reason- 
ing. This destroys the monotonic convergence property of the Krotov algorithm. However, 
useful results can be obtained by storing in memory during the process the best field ob- 
tained so far. The field that is stored after a predefined number of iterations is the result 
of the optimization process. 

2.4.2 The method of Degani et al. : A non-diagonal penalty term 

Degani et al. [12] developed a method of restricting the spectrum by a more complicated 
penalty term, J penal- This enables the inclusion of more information on the desired field. ^ 
As an introduction, we first mention a simpler use of the penalty factor, in order to 
control the properties of the field (see, for example, [2]). A time dependent penalty factor, 
a{t) can be introduced, in order to control the time-shape of the field. The modified J penal 
is given by: 



^The formulation of the control problem in is discrete, due to the introduction of control by a 
piecewise constant field, in the same paper. To avoid complications and inhomogcneity in the present 
text, we reformulate the problem in a continuous context. However, note that there is indeed a need for a 
discrete formulation, when dealing with the numerical method of maximization. 




(2.46) 




(2.47) 
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For instance, we can force a Gaussian envelope on the field profile, by choosing an appro- 
priate a(t). 

The a{t) function is one dimensional, and contains only one dimensional information. 
Degani et al. use a two dimensional penalty term, which contains greater wealth of infor- 
mation about the properties of the desired field. We introduce the two dimensional penalty 
function f3{t,t')] the new penalty term takes the form: 

JpenaiW)] = - f [ e{t)l3{t,t')e{t') dt dt' (2.48) 

An alternative formulation is in the operator language, where B stands for an arbitrary 
linear operator: 

JpenaiW)] = - I e{t)Be{t) dt (2.49) 
Jo 

Be(t) = [ /3{t,t')e{t')dt' (2.50) 
Jo 



In this formulation, the role of I3{t,t') is more apparent. 

The B operator (or matrix, in the original context; see footnote 2) used by the authors 
is the following: 

B = OigoodPgood + OibadPbad (2.51) 
< abad > 



Pgood is the projection operator on the subspace of desired control functions; Pbad is the 
projection operator on the subspace of undesired control function, where: 

Pbad I Pgood 

(I is the identity operator). The a's have a role similar to a penalty factor. The idea is to 
"penalize" the optimization process for undesired control functions and to encourage the 
appearance of desired control functions. In this case, the desired functions are fields in the 
desired frequency domain. In this way, we can control the spectral properties of the field. 
(Another suggestion for B is introduced by the authors; we will not discuss it here.) 

The resulting Euler-Lagrange equations are the same as in Sec. 2.1, apart from the 
replacement of Eq. (2.21) by the following: 

Be(t) = / /3(t, t')e{t') dt' = -Im(x(t) m) (2-52) 

The new equation does not have the special property mentioned in subsection 2.3.3 - 
e{t) does not depend only on other objects at the same time t, but also in the field at 



•^The authors include an additional term in B: al, maybe for convenience. We didn't include it, because: 
I = Pgood + Pbad, SO this term is unnecessary. 
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all other time points. This prevents the possibility of using the Krotov method in the 
way it appears in subsection 2.3.3. Instead, the authors present a Krotov-like method, 
which requires the use of a discretized version of Eq. (2.52) (of course, discretization of 
the problem is required anyway for a numerical solution). Let us define a time grid (the 
authors use an equidistant grid; here we define a grid which is not necessarily equidistant): 

ti,t2, ■ ■ ■ ,tN 

Now, we approximate Eq. (2.52) by a discretized form of the equation: 

N 

TY,Pitut,)e{tj)w, = -Im(x(t.) |/i|^/'(t^)) (2.53) 



where Wj is the integration weight for time point tj. Let us define: 
Rearrangement of Eq. (2.53) gives: 



hj = Twj 



, . _ -Im(x(t^) lAl jjitj)) - Ei=i/^(^i^^j)e(^i)^j - EjLm/^(^"^i)^fe-)^i ,^ . 

This equation is the basis for the update rule of e{ti) in the Krotov-like algorithm. The 
idea is to compute in each time point the field with the newest values available of the 
variables. Eq. (2.44) is replaced by: 



I3{ti,ti)hi 

(2.55) 



Eq. (2.45) is replaced by: 

-Im(x('=)(t.) V^('=+^)(t.)> - j:T=\m,h)e('''-'Kh)h, - Ef=m/5(^-^.)e~^'^/^, 



iu) 



(3{ti,ti)hi 

(2.56) 



Unlike the standard Krotov algorithm, this algorithm lacks the self consistency of the 
sequence of e(t), because of the use of data from another sequence to compute e(ti). Nev- 
ertheless, the authors report monotonic convergence of the algorithm. However, a conver- 
gence analysis is not supplied. 



2.4.3 The method of Skinner and Gershenzon: controlUng a Ust 
of frequency terms 

Skinner and Gershenzon [13] present an approach of controlling the properties of the field 
by defining the field as an analytic function of a desired form with several adjustable 
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parameters. The parameters are the optimized variables, instead of e{t). In our case, the 
field is defined by a cosine series: 

N 

e{t) = ^ a„ cos(n Aw t) (2.57) 

n=0 

where Aw is defined by the resolution of w, at a time interval of T: Aw = tt/T. The 
maximal frequency that may be present in the resulting field is: A^Aw. The a„ parameters 
are optimized to maximize J. 

The optimization of different variables requires an alteration in the maximization con- 
dition; Eq. (2.14) is replaced by the following set of equations: 

dJ 

= n = 0,l,...,Ar (2.58) 

dan 

The maximization problem that is dealt with by the authors is somewhat different from 
ours; hence, we avoid presenting the resulting Euler-Lagrange equations. 

The optimized variables are not defined in the time domain, like e{t). This rules out 
the possibility of using the Krotov method. The authors use a first order gradient method. 



2.5 OCT of a time dependent dipole moment 

In order to generate an emitted radiation of a desired frequency, we have to be able to 
control the frequency of the oscillations of the time dependent dipole moment expectation 
value: 

{m) = m)mm) 

Hence, the development of methods of controlling {fL){t), are of utmost importance in the 
task of harmonic generation control. 

The main problem in controlling the path of {fi){t) by OCT is that this problem is 
not a maximization problem, when formulated in the terms of {fi)(t). Hence, we obviously 
cannot use the method of Sec. 2.2 by simply setting: 0(t) = fl. 

Serban, Werschnik and Gross [11] propose to solve this problem by maximizing an 
alternative time dependent operator, instead of dealing with fi directly. Then, the method 
of Sec. 2.2 is applicable. Consider a case of a quantum system with one spatial variable, x, 
and a charge: q = 1 for convenience. In this case: /i = X. The time-dependent operator 
is: 

6{t) = 5 [x - ^{t)i\ (2.59) 
where S{x) is the Dirac delta function. Practically, it is approximated by a sharp Gaussian: 



6{x) ~ \j — exp{—bx 
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where b is large, ^{t) represents the desired path of {fi){t). When (|(5[X — is 

maximized, the density of the wave function in the neighbourhood of x = ^(t) is also 
maximized. Hence: 

(/i)(t) = (x)(t)^e(t) 

By specifying the desired path ^{t), we will hopefully be able to control the path of {fi)it). 

In this approach the problem is treated semi- classically by considering the wave- 
function as a localized object. This restricts considerably the possible mechanisms for 
varying along a predefined path. Moreover, nothing ensures that the semi-classical 

mechanism is at all possible for a given path, while other mechanisms may be possible. 

The method was employed rather successfully for a path that was known in advance to 
be possible. 

If the method is to be used for harmonic generation, ^(t) has to be defined as an 
oscillating function with the desired frequency. This has to be combined with the restriction 
of the forcing field spectrum. 

The authors reported (in 2005) on a research in high harmonic generation control. We 
assume, that the brute force method (introduced by two of the authors, Werschnik and 
Gross) was employed for the restriction of the forcing field spectrum. The results of this 
research were not published. We know from private communications that the reason is 
that this attempt has failed. 

A possible explanation for the failure of this attempt is the one mentioned above: the 
method enforces a semi-classical mechanism, which is not necessarily possible. If this is the 
reason for the failure of this method, another formulation for controlling (fi) (t) is necessary 
for the harmonic generation problem. Another possible explanation is that the brute force 
method is applicable only for simpler targets, like in Sec. 2.1. If this is the reason for the 
failure of this attempt, other methods for restricting the field spectrum may be attempted, 
combined with the formulation described in this section. 
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Chapter 3 

The new method 



In this chapter, we present the new method for harmonic generation by QOCT. We present 
a new formulation for the maximization problem, and discuss the numerical methods that 
can be applied in this formulation. 

As an introduction, we point out the main difficulty in formulating the harmonic gen- 
eration problem in the context of QOCT. The quantum dynamics is formulated in the 
time domain. In the common QOCT formulation the quantum dynamics is represented 
by the time-dependent Schrodinger equation. In the regular control problems (that were 
presented in the previous chapter. Sec. 2.1, 2.2), the maximization requirements are related 
to one, or many time points, each with its own requirement. Hence, the whole problem 
is naturally formulated in the time domain. As for the harmonic generation problem, the 
requirements on the forcing and emitted fields are related to the frequency domain. Such re- 
quirements are naturally formulated in the frequency domain, but not in the time domain. 
A requirement on a single frequency is related to the whole time sequence simultaneously 
and not to a single time point. When we deal with requirements on a spectrum sequence, 
with many frequencies, the situation is even more complex. 

The difficulty is to create a unified formulation that takes into account the quantum 
dynamics and satisfies also the frequency requirements. We have already seen several ways 
to deal with this difficulty, in the discussion on the methods for restricting the forcing field 
spectrum (Sec. 2.4). 

The approach that was adopted here to deal with this problem is very simple: Each part 
of the problem is handled in its natural domain — the quantum dynamics is formulated 
in the time domain, by the time-dependent Schrodinger equation, while the frequency 
requirements are formulated in the frequency domain, by using appropriate functionals. 
The switch between the time and frequency domains is made by the cosine transform (we 
adopted the cosine transform, instead of the more commonly used Fourier transform, for 
reasons that will be mentioned). This makes the resulting equations to look somewhat 
cumbersome; nevertheless, it is possible to recognize the meaning of the resulting forms. 

The treatment of the problem of harmonic generation is divided in this chapter into the 
two parts of our task, mentioned at the beginning of the previous chapter: in Sec. 3.1, we 
propose a new formulation for handling the general problem of imposing a restriction on 
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the forcing field spectrum. Sec. 3.2 deals with the numerical methods that are appropriate 
for this formulation. In Sec. 3.3, we introduce a new maximization functional for the 
requested spectrum of the dipole, and present the full maximization problem for harmonic 
generation. 

3.1 New formulation for imposing restrictions on the 
forcing field spectrum in QOCT 

In this section, we present a new formulation for dealing with the general problem of 
restricting the spectrum of the forcing field. In Sec. 3.3, this formulation will be used for 
the harmonic generation problem. 

As will be seen, there is a close relationship between all the methods of Sec. 2.4 and 
this formulation. 



3.1.1 Frequency dependent penalty factor 

Before we start, we introduce the cosine transform, that will be used frequently throughout 
this chapter. We denote it by the symbol C. The cosine transform of an arbitrary function 
g{t) is defined as: 

C[g{t)] = J- [ g{t)cos{ujt)dt (3.1) 



The function g{uj) is defined as the result of the operation of the cosine transform on g(t): 

-g{u) = C[git)] (3.2) 

We call it: "the cosine transform of g(t)" . 

The inverse of the cosine transform is the inverse cosine transform, that will be denoted 
as C-^: 

C-'[g{uj)] = ^J^ g{uj)cos{ujt)du = g{t) (3.3) 

Note that the inverse cosine transform is the same as the cosine transform, except the 
variable names. 

We prefer the cosine transform as a spectral tool over the more commonly used Fourier 
transform for two reasons: 

1. The cosine transform of a real function is also a real function. This is not true for 
the Fourier transform of a real function. By using the cosine transform for a real 
function, we avoid complications due to the appearance of a complex transformed 
function (we have already encountered such complications in the previous chapter. 
Sec. 2.1, with l^pit))). 
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2. When handling with the resulting equations numerically, the cosine transform is 
approximated by a discrete cosine transform (DCT). In our case, the DCT is prefer- 
able over the discrete Fourier transform (DFT), because the "ringing" phenomenon 
(known also as "aliasing") is much less pronounced in the DCT; hence, we can avoid 
the use of "windows" ^ 

Now we return to our problem. The object of main interest for us is e(ci;), which 
represents the spectrum of the field. We want to find a way to affect its general shape. 

Note that the integration domain of the cosine transform does not have to be infinite 
in this case. The reason is that e{t) is defined only in the time interval: < t < T; hence, 
e{u) can be written finite time integral: 

2 '•^ 



e{Lo) = y- J cos{ujt) dt (3.4) 

The same is true for the integration over a;, in the inverse cosine transform: In practice, 
e{uj) is negligible over some value of cj. Moreover, when treating the problem numerically 
the maximal frequency possible is limited by the resolution of the time points by the 
Nyquist- Shannon theorem; according to this theorem, the maximal frequency cannot ex- 
ceed vr/At, where At is the distance between neighbouring time points. Anyway, we can 
define the maximal relevant frequency in the system as f2. The inverse cosine transform 
can be written with VL as the upper limit, instead of oo: 

e{t) = \ — / e{(jj) cos{(jjt) du (3.5) 
V Jo 

In order to restrict the spectrum of the forcing field, we introduce a new penalty func- 
tional, Jpenah which is formulatcd in the frequency domain instead of the time domain. We 
use a frequency dependent penalty factor: 

JpenailK'^)] = " / Ci{uj)f {u) duj a{uj) > (3.6) 

Jo 

The idea is to choose very large values for a{uj) for the undesirable u values, and small 
values for the desirable values. Later in this section, a{oj) will attain a more precise 
meaning. 

It is more convenient to replace the optimized function e{t) by e{uj). Accordingly, we 
replace the condition for an extremal in Eq. (2.14) by the following one: 

(3.7) 



^The ringing results from discontinuities in the extended periodic function, represented by the truncated 
discrete Fourier series, or in its derivatives. This topic wiU be discussed further in Subsection 3.3.3. See 
also [20], for a more detailed explanation on the advantage of the DCT over the DFT. 
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When taking the functional derivative, we note that only Jpenai and Jcon are dependent 
on e{u) (when treated as independent of I'lpit)) and \x{t)): according to the Lagrange- 
multiplier method). This makes the present treatment suitable for diverse kinds of problems 
(the same is true for the methods mentioned in the previous chapter). The dependence of 
Jpenai IS simple. The dependence of Jcon is through e{t) in the time-dependent Hamiltonian: 

H(t) = Ho-/le(t) = Ho-/t ^y^^ e{uj) cos{ut) du}j (3.8) 

The resulting Euler-Lagrange equation for e{uj) is (the full derivation is given in App. A): 

[2 j^lm{xit) \il\m)cosiut)dt 

e{ij) = -A / — y — (3.9) 

V TT a[UJ) 

This equation replaces Eq. (2.19) for e(t) in regular control problems. In our problem, e{t) 
can be easily computed from e{u), by the inverse cosine transform. 
The meaning of Eq. (3.9) will be more apparent if we define: 

a{u;) = -f^ (3.10) 
a>Q (3.11) 



n 

f,{uj)du = l (3.12) 

a is a positive constant. /e(u;) contains the dependence of a{u) on u. This function is 
always positive, due to the conditions in Eqs. (3.6), (3.11). The normalization condition 
in (3.12) makes /e(w) and a well defined. 

We also define, for convenience, the following function of t: 

a 

Note that the RHS is just the expression for e(t) in Eq. (2.21), for the regular control 
problems, where a from (2.21) is replaced by a. 
Now, let us write Eq. (3.9) as follows: 

2 '■^ 



e(u) =fJu)\ — / rjit) cos(ut) dt 
V vr Jo 

=aco)C[v{t)] (3.14) 
The expression for e(t) in our problem can be written as: 

e{t) = C~' {Uu;)Cm]} (3.15) 
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Now, the meaning of the resulting Euler-Lagrange equation is obvious: We start from the 
field in the regular control problems, represented by T]{t); a has the role of a global, constant 
penalty factor. We transform this field to the frequency domain. Then, we multiply it by 
/e(a;), which has the meaning of a filter function. The undesirable frequency components 
are filtered out, and an envelope function can be forced on the profile of the spectrum. 
Then, we transform the resulting spectral function back to the time domain, and get the 
e{t) sequence. 

The simplest choice for /^(a;) is the rectangular function. If the allowed frequencies are 
in the interval: [wmm, ^max], then the rectangular function is: 

<UJ < Umin 

fe{^) = { 7, ^rnin <0J <UJ. 



UJmnx < (jJ 



max 



-U{UJ - UJ„iin)u{uj.max - Uj) (3.16) 



^max 

where u{x) is the Heaviside step function, defined as: 

uix) = < 3.17 
^ ^ [1 0<x ^ ^ 

With this choice for fe{u}), we get a complete filtration outside the interval. To be more 
strict, the value of /e(cu) cannot be absolutely — otherwise, a{uj) is undefined. It is more 
precise to refer to the limit, in which outside the interval fei^j) tends to 0, and then, a{u) 
tends to oo. Nevertheless, Eq. (3.16) can be used as is for practical purposes, together 
with Eq. (3.14). 

Other choices for fe{oo) are possible; for instance, we can choose a Gaussian function, 
or a "hat function" (see Fig. 3.1, for an example of a hat function). These functions can be 
used when we want a "smooth" filtration. We can also choose /e(w) with a special shape 
to enforce a desired envelope shape on the profile of e(w). 

The division of the quantity l/a{u) into the fraction of two quantities: a and feijjj), is 
instructive conceptually; however, it has no meaning for practical use. The reason is that 
the exact value of a has no precise meaning — it has to be determined by a trial and error 
process, and cannot be known in advance. Hence, there is no advantage in normalizing 
/e(cj). In practice, it is convenient to use the function /e(w) instead, defined as: 

/eM = ^ (3.18) 

For practical use, Eq. (3.9) is written conveniently as: 

-e{u) = Ucu)C [-lm{x{t) |A| m)] (3-19) 



27 




0, 







0.5 



1.5 



2 



X 



Figure 3.1: An example of a "hat function" — the function: y{x) = sech[20(a; — 1)^] 

3.1.2 The relation between the new formulation and the existing 
methods 

Now we show the relation between the new formulation and the existing methods for 
restricting the spectrum of the forcing field. The formulation of Degani et al. is shown 
to be equivalent to the new formulation; that of Skinner and Gershenzon is shown to be 
a special case of it. There is also a close relation between the new formulation and the 
method employed by Werschnik and Gross. 

We start with the method of Degani et al. . In our treatment of the new Jpenau we pre- 
ferred to work in the frequency domain, with e{uj) as a variable; in this way, the frequency 
requirements are expressed more naturally. However, it is also possible to translate the 
whole problem to the time domain, with e{t) as a variable. We express e{u:) in Jpenai as a 
cosine transform of e{t) (Eq. (3.4)). The resulting Jpenai is expressed in the terms of e(t): 



Jpenai[e{t)] = - a{uj)\ - e{t) cos{ujt) dt \ - e{t') cos{ujt') dt' doo (3.20) 





After rearrangement, Eq. (3.20) can be written as follows: 




(3.21) 




n 



(3.22) 



Note that Eq. (3.21) is the same as (2.48) in the formulation of Degani et al. . 
Alternatively, Eqs. (3.21), (3.22) can be written as: 




(3.23) 



(3.24) 
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Eq. (3.23) is the same as (2.49). The meaning of the operator B in Eq. (3.24) is apparent: 
it transforms e(t) to the frequency domain; then each frequency component is multiphed 
by its own penalty factor, a{uj); the resulting function of u is transformed back to the time 
domain. 

The operator B proposed by Degani et al. (Eq. (2.51)) can be written in the form of the 
operator in Eq. (3.24), with an appropriate a{uj). Suppose that the "good" field functions 
are in the frequency domain: [0, Umax]'-, then, the projection operators from (2.51) are 
defined by the following expressions: 

Pgoode(^) = {u{uj.max " uj)C[e{t)]} 

Pbade{t) = {U{UJ - UJmax)C[e{t)]] 

(Of course, another spectral transform can be chosen instead of the cosine transform.) The 
operator B from Eq. (2.51) can be written as in (3.24), with the following a{(jj) (using the 
linearity property of the cosine transform) : 

a{uj) =agoodU{uJmax - Uj) + ahadU{u} - Umax) 
agood <UJ <U. 



■'max 



(3.25) 

dbad ^max < ^ 



Note that Degani et al. use an extended definition of a(cj), where it can attain negative 
values. Then, the interpretation of a{uj) that was given here is not appropriate. We have 
decided to avoid using negative values to prevent the possibility of singularities in the 
hypersurface of J (see Sec. 2.1 in the previous chapter). 

We see that the formulation of Degani et al. is equivalent to ours. The difference 
between the two approaches is that we work in a basis in which B is diagonal — the 
frequency basis, while Degani et al. work in a basis in which it is non-diagonal. It follows 
that our treatment is one-dimensional, while the treatment in the formulation of Degani 
et al. is two-dimensional. This is due to the fact that the requirements on the frequency 
are naturally formulated in the frequency domain, and the resulting expression for Jpenai 
is considerably simpler. The resulting Euler-Lagrange equation for the field is also more 
naturally expressed in the frequency domain — compare the integral equation (2.52) for 
e(t), with Eq. (3.19) for e{u). 

One important disadvantage of the approach of Degani et al. is that a complete fil- 
tration, like the one achieved by the rectangular fe{uj), is impossible. The reason is that 
a{u) diverges, and there is no way to create the corresponding /3(t, t'). The problem is not 
restricted only to a rectangular filter function — any function that decays very rapidly (as 
in the case of an exponential decay, e. g. in a Gaussian function, or in the hat function of 
Fig. 3.1) is not suitable, because a{u) attains very large values that are not acceptable 
numerically. Another problem with rapidly decaying functions is the very sharp shape of 
the resulting <y{uj), which is not suitable for numerical treatment. For instance, in practice, 
the integral in Eq. (3.22) cannot be performed numerically, if a{uj) is too sharp. In Sec. 3.2 
we will discuss another numerical problem which has a similar origin. 
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The two formulations result also in different numerical approaches for treating the 
problem; this will be discussed in Sec. 3.2. 

Now we show that the formulation of Skinner and Gershenzon is a special case of the 
new formulation. 

Consider the following filter function: 



N 



2N 



— 5{uj — n Aw) 



(3.26) 



n=0 



The factor 2/{2N + 1) is a normalization constant (note that integration on S{uj) for the 
n = term gives only 1/2, because the integration domain in Eq. (3.12) starts from 0). 
Let us derive the expression for e(t), for this filter function. From Eq. (3.14) we have: 



e{uj) = f,{uj)ri{uj) 
Using Eq. (3.5), we get the following expression for e(t): 



(3.27) 
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Comparing Eq. (3.28) with Eq. (2.57), we can recognize that we got the same expression 
for e{t) as that was used by Skinner and Gershenzon, with: 
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Eq. (3.29) is the Euler-Lagrange equation for the a„'s. We see that the problem of op- 
timizing a continuous function was reduced, with this choice of /e(w), to the problem of 
optimizing a discrete set of parameters. 

It is important to note, that in numerical computation the general continuous problem 
of the new formulation is approximated by a discretized version of the problem. This 
amounts to replacing the continuous cosine transform by a discrete cosine transform, which 
is a cosine series. If we also use a rectangular /^(w), the problem solved numerically is 
essentially the same as in the formulation of Skinner and Gershenzon (see App. B). The 
present formulation is still more general, because of the possibility of using other forms 
of /e(ci;), for a smooth filtration, or in order to affect the shape of the spectrum. These 
advantages are applicable also for a discrete version of the problem; starting from the 
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general continuous problem of the new formulation, we may use a generalized version of 



where the fe are a set of predefined constants which determine the desired shape of the 
spectrum in a similar manner to the continuous case. The derivation of the Euler-Lagrange 
equations is the same as above. 

The brute force method of Werschnik and Gross is not equivalent to the new formula- 
tion. However, the new formulation provides some justification to this method. 

Compare the field that is used to propagate \ip(t)) in the brute force method (Eq. (2.46), 
together with Eq. (2.44)), to the field from Eq. (3.15). The two equations are very similar. 
The only differences are: 

1. The Fourier transform from Eq. (2.46) is replaced by a cosine transform in Eq. (3.15). 

2. a of the regular formulation is replaced by a. 

3. f{uj) from Eq. (3.15) is replaced by f\{uj). 

According to the interpretation given here for a and fe{u}), the meaning of the two equations 
is essentially the same. 

However, the field used to propagate \x(t)) in the brute force method (Eq. (2.44)), is 
not justified by the new formulation — the function r]{t) does not have the interpretation 
of the field appropriate for this problem, before it is filtered, as in Eq. (3.15). 

3.2 Numerical methods for the new formulation 

In this section, we discuss a few numerical methods that can be suggested for maximizing 
J, with the Jpenai presented in Sec. 3.1. 

Usually, the Krotov method (or one of its variants) is considered the preferable method 
for solving QOCT problems. This is because of its monotonic convergence property and 
the relative ease of its application. However, the Krotov method is not applicable for 
the new formulation; the reason is, that e{t) from Eq. (3.15) lacks the property of being 
dependent only on functions of the time t. This makes it necessary to look for other 
numerical methods, that are applicable and efficient when employed for this problem. 

There are two possible approaches for numerical solution to our problem: 

1. We can use the Krotov- like approach, presented in the previous chapter (Subsec- 
tion 2.4.2). In this approach, we utilize the fact that the sequence of \ip{t)) or |x(^)) 
does not have to be determined simultaneously in all times, for updating the field 
during the propagation. For applying the Krotov-like approach it is necessary to 
formulate the problem in the time domain, in the way that was presented in the 
previous chapter. Adopting this approach means, essentially, turning to the existing 



Eq. (3.26): 
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method presented by Degani et al. , with a somewhat different approach for the B 
operator. 

2. We can use the more natural formulation of the problem in the frequency domain, 
and employ general methods for nonlinear optimization problems (the time-domain 
formulation is not preferable for this approach, because the resulting update rules 
are more complicated, and more expansive numerically). 

Subsection 3.2.1 deals with the first approach; Subsections 3.2.2-3.2.3 represent the second. 
3.2.1 The Krotov-like method of Degani et al. 

When we tried to implement the Krotov-like algorithm, we encountered problems. The 
algorithm has been employed for the harmonic generation problem, with Jmax that has not 
been introduced yet, and will be introduced in Sec. 3.3. Anyway, the problems that we 
have encountered with seem to be general problems with this algorithm. Other groups also 
report (in private communications) on similar difficulties when applying the algorithm to 
the common control problem of Sec. 2.1. 

The algorithm was found to be causing a numerical instability during the propagation 
process. The origin of the instability was found to be the fact that /3(t, t') is very oscillatory 
in nature. The amplitudes of the oscillations become larger, when a{u:) is chosen to 
attain larger values for undesirable frequency components. Larger oscillations increase the 
numerical instability. 

It is necessary to mention here that the harmonic generation problem is a difficult 
problem. One reason is the typically small amplitudes of the higher frequency oscillations 
of the dipole that can be achieved by the low frequency field. Another reason is the difficulty 
in achieving a maximum, probably because of a complex structure of the hypersurface of 
the optimization functional J. It follows, that there is a need for very large values of ot{u:), 
to achieve a satisfactory filtration; otherwise, the optimization algorithm will prefer an 
easier path for increasing J, using high frequency fields. The values of a^ad that were used 
by Degani et al. in [12] are several orders of magnitude smaller than the values that are 
required for the harmonic generation problem. When using a{(jj) values of the same order 
of magnitude that was used in [12] the numerical difficulty was found to be acceptable; 
however, such values do not satisfy the requirements of the harmonic generation problem. 

The origin of this general problem with the algorithm can be considered to be the 
unnatural formulation of the problem in the time domain. This formulation makes it 
necessary to use the time-dependent, two-dimensional penalty function f3(t,t'), instead of 
the much simpler a{uj). The relatively simple frequency requirements of the problem are 
expressed in the time domain only by a complicated, non-smooth structure of P{t, t'), with 
large oscillations of very high frequencies. 

Another problem that we encountered (even with relatively small values of (y{u})) is 
related to the convergence characteristics of the algorithm: the algorithm has shown mono- 
tonic convergence characteristics only for forcing fields with relatively low intensities. When 
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the forcing field intensity became larger, the monotonic convergence characteristics were 
destroyed. We have not tried to find the origin of this problem. Possibly, the problem 
is just a result of inaccuracies in the propagation process due to the numerical instability 
mentioned above. On the other hand, the fact that this problem was reported also by 
other groups, together with the absence of available convergence analysis, allows to raise 
the suspicion that this Krotov-like algorithm is not always monotonically convergent (al- 
though it may show frequently monotonic convergence properties, as shown in [12], a fact 
that may be useful). 

It is noteworthy that we tried to use another Krotov-like algorithm, based on updating 
the field in every time step of the propagation by Eq. (3.15). This update rule is even less 
self-consistent than that used by Degani et al. . Although this algorithm shows sometimes 
monotonic convergence characteristics, usually it is even worse than the algorithm of the 
naive approach (Subsection 2.3.1). 

The failure of the attempt to apply a Krotov-like algorithm for our problem, makes it 
necessary to turn to the second approach mentioned above. 



3.2.2 The BFGS method 

A natural choice for a general optimization algorithm is a quasi-Newton algorithm, because 
of the good convergence characteristics near the maximum. In this subsection, we deal with 
the BFGS method, which is the most commonly used quasi-Newton method. 

In the new formulation for imposing restrictions on the field spectrum, the update 
rule is more conveniently formulated for e{u:). Hence, we have to rewrite all the relevant 
equations for the quasi-Newton methods (Subsection 2.3.2) with e(a;) as the variable. The 
gradient with respect to e(a;) is^: 

J = = -2 {a{uj)-e{uj) + C [Im(x(t) ^(t))]} (3.31) 



The Hessian is defined as: 



5V 

s{uj,uj') = (3.32) 
de{u)de{uj') 



^Eq. (3.31) raises a problem: when using a rectangular /e(w) (Eq. (3.16)), a{uj) diverges for the unde- 
sired ui values. However, for these values, eiuj) is certainly 0, so we can set in advance: 

e(a;) = 0, for: UJ < LOmin, > i^max 

without including these components of the field in the optimization process. This is equivalent to computing 
e(t) only from frequency components in the domain: [wmm, ^^max], by changing the integration domain in 
Eq. (3.5) to this domain. The integration domain in (3.6) has to be changed accordingly. A similar problem 
exists for rapidly decaying /e(i^) functions, as was mentioned in another context in Subsection 3.1.2. The 
solution for this problem is similar: e(cj) is set to when the corresponding a(cj) is large enough, so e(a;) 
can be assumed to be negligible (note that such a solution is impossible when dealing with the time domain 
formulation). The same problem exists for the computation of J during the line-search, and the solution 
is similar. 
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The operation of the Hessian on an arbitrary function of cu, g{u), is defined by: 

Sg{u)=[ s{uj,uj')g{uj')duj' (3.33) 
Jo 

The update rule in Eq. (2.43) is replaced by the following equations: 

e'"-\cu) = e^'^cu) - K S%^-'[V,(^^)jf^ (3.34) 
e*"'^'(t) =C-^[e*""V)] (3.35) 

After the line-search was performed, and an optimal K was found, e{u) is updated; we add 
to Eq. (2.34) the additional update: 

e('=+i)(c^) = (a;) (3.36) 

The BFGS method was implemented by using the f minunc function from the "Opti- 
mization Toolbox" of MATLAB (see [5] for details). The method was applied successfully 
to simple problems of a two-level-system (TLS), with the Jmax of sections 2.1 and 3.3. 
Unfortunately, the method was found to be much time consuming in this case, so this 
option is not convenient for larger-scale problems. This raises the motivation to find an 
alternative method. 



3.2.3 The relaxation method 

The "relaxation method" is a general method for "helping" an iterative process to converge. 
To our knowledge, it has not been employed yet for QOCT problems. In this subsection, we 
start from the description of the relaxation method, in the context of the simpler problems 
of the previous chapter, sections 2.1, 2.2. A mathematical justification to the method is 
given. Numerical results for the simpler problems are presented . Then, the method is 
discussed in the context of the new formulation of Sec. 3.1. Numerical examples will be 
given for problems with the Jmax of Eq. (2.1). 

The relaxation method, in the context of regular QOCT problems, is based on the 
following update rule: 

^new^^^ = Ke^^it) + (1 - K)e°^^{t) 0<K <1 (3.37) 



^ELr 

. _ lm{x{t)\jl\m) 
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(3.38) 



e (t) is the field from the Euler-Lagrange equation (2.21), using the |V'(^)) and |x(^)) 
sequences, propagated by e°''^(t). Actually, this is the new field, according to the update 
rule in the naive approach (see Subsection 2.3.1). The idea is to "mix" the old solution with 
the new solution of the naive approach. is a parameter, which determines the weights of 
the old and new solutions. Its value has to be determined by a trial and error process. Note 
that if we take: = 1, we return to the naive approach update rule. Although the naive 
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approach seldom converges, the relaxation method may be convergent, by an appropriate 
choice of K. Usually, the value of K should be decreased during the optimization process. 

We propose the following scheme for the implementation of the relaxation method, for 
the standard QOCT problems: 

1. Guess a field sequence: e*^°^(t). 

2. Guess an initial value for K. 

3. Propagate \il)^'^\t)'^ forward from t = to t = T, with e''°^(t). 

4. Calculate with \ip^^\t)) and e^'^\t). 

5. (k = 0) 

6. Repeat the following steps, until convergence: 

(a) Set \x^^\T)'^ according to the suitable boundary condition, using |^/^*^'^)(T)). 

(b) Propagate \x^^\t)) backward from t = T to t = 0, with e^''\t). 

(c) Do the following steps, and repeat while J**"*"^' < J^'^); 

i. Set a new field, using Eq. (3.37): 



e*"'*'(t) = K 



Im( x^''Kt) 



a 



+ (l-ir)eW(t) 



(3.39) 



\t)) and e*"'''(t). 



ii. Propagate |'?/;*"'^'(t)) forward from t = to t = T, with e*""'(t). 

iii. Calculate J^""^ with |^*"«'^+^\ 

iv. If J*""' < J^'^), then set: K = K/2 
(d) Update all the variables: 



trial 



J- 



trial 



(3.40) 



(k = k + 1) 



During this procedure, the value of K decays exponentially with the number of its 
updates. Most frequently, it decays in the former iterations, and remains constant during 
the rest of the process. 

We denote the initial guess for K in step 2 as: Ki. Usually, i^j = 1 is a good choice. 
However, often it may be advantageous to choose a smaller value, for the following reasons: 

1. When a is very small, a numerical instability in the propagation process might ap- 
pear, unless K is small enough (this numerical problem exists also in the Krotov 
method). 
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2. Large K may be advantageous in the former iteration; however, after many iterations, 
we may get a better result with a smaller Ki. 



Now we are going to give a mathematical justification for the relaxation method. We 
show, that the relaxation method is actually a quasi- Newton method. As we mentioned 
in the previous chapter (Subsection 2.3.2), in practice, the Hessian (Eq. (2.38)) cannot 
be computed in every iteration, because of the complex dependence of \ilj{t)) and |x(^)) 
on e(t) (which are present in the expression for the gradient, Eq. (2.31)). Now we claim, 
that we can get a useful approximation to the Hessian by ignoring this dependence, and 
treating \ip{t)) and |x(^)) independent of e{t). This is equivalent to the approximation 
of the Hessian by the Hessian of Jpenai- To be more general, we include the possibility of 
a time-dependent penalty factor, a{t) (Eq. (2.47)). The approximation for the Hessian is 
defined by the following equations: 



S{t,t') ^ Sap{t,t') 



52 J, 



penal 



-2a{t')5{t' - 1) 



6e{t)6e{t') 

The operation of the Hessian on an arbitrary function of t, g(t), is defined by: 

^apg{t) = -2 f a{t')5{t' - t)g{t') dt' = -2a{t)g{t) 
Jo 



(3.41) 
(3.42) 



(3.43) 



We see that S^p is diagonal in the time basis. Now, it is a trivial matter to get the inverse 
of Sap-. 

Sap'ait) = (3-44) 

Let us write the expression to the update rule in the quasi-Newton method (Eq. (2.42)), 
with the new Sap, using (2.31) and (3.44): 



e''"*(t) - K- 



2a{t) 

--Ke^'^it) + (1 - K)e°'\t) 



a{t)e"'\t) + lu.{x(t) V^m)) 



e(t)=eoW(4) 



which is the same as (3.37). Note, however, that K has a slightly different definition in a 
quasi-Newton method, because it can attain, in principle, values greater than 1. 

In the regular Newton method, with the exact Hessian, K is absent, which is equivalent 
to setting: K = 1 (see Eq. (2.40)). The fact that the naive approach, with K = 1, 
is unsuccessful, indicates that the approximation of Eq. (3.41) is not a very good one. 
Anyway, this approximation proves itself to be useful for determining the direction of 
search, in the context of a quasi-Newton method. 
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Note that when a is a constant, the operation of 



on the gradient vector is just a 



multiphcation by a negative constant. The resulting direction of search is the direction of 
the gradient, as in the first order gradient method (see Eq. (2.32)). However, the gradient 
vector, which determines the direction of search, is multiphed here by the constant l/(2a). 
This may have an effect on the numerical search procedure. 

This new interpretation of Eq. (3.37) raises the possibility of using the common scheme 
for the quasi-Newton methods, with a line-search in every iteration, instead of the decay 
of K proposed above. To distinguish between the two possibilities, we will call the scheme 
that was proposed above: "the decaying K scheme" , and the one that employs a line-search: 
"the line-search scheme" . 

First, we discuss the numerical results of the decaying K scheme. 

To demonstrate the efficiency of the relaxation method, we compare the convergence 
curve of this method with that of the Krotov method. In the convergence curve, we 
plot J vs. a variable which is a measure of the numerical effort; here it is the number of 
propagations. We have chosen two simple state-to-state problems, with Jmax of the form 
of Eq. (2.1), where the projection operator from Eq. (2.3) stands for O. 

The first problem is a TLS problem. The problem is to find a field for a transition from 
the initial state — the ground-state, to the target state — the excited state. In our TLS 
examples, we take the dipole moment operator to be the x Pauli matrix: 



The rest of the details of this problem, and the resulting curve, are presented in Fig. 3.2. 

The system of the second problem is a one- dimensional harmonic oscillator. The initial 
state is the ground-state. The target state is a "shifted" ground-state. The details and the 
resulting curve are presented in Fig. 3.3. 

The results show, that in these simple cases, the relaxation method is definitely ad- 
vantageous over the Krotov method. An additional research is required to determine the 
efficiency of the relaxation method, compared with the Krotov method, in more complex 
problems. This is beyond the scope of our research. 

We have also tried the line-search scheme. It has been implemented using f minunc. We 
found that this method is much less expansive than the BFGS method, and it converges to 
better solutions. However, it is considerably more expansive than the decaying K scheme. 
Moreover, somewhat surprisingly, the resulting solutions, after convergence, are of smaller 
J values, compared with the decaying K scheme solutions. Hence, the decaying K scheme 
was adopted in the present work, and is used in all the examples throughout this text. 

It is noteworthy, that after the decaying K scheme has converged to some solution, it 
is possible to get some improvement using the line-search scheme, with this solution as a 
first guess. We observed also, that a further improvement can be achieved, by using the 
BFGS method after this stage. 

To employ the relaxation method to the new formulation of restricting the forcing field 
spectrum we have to reformulate the update rule in the frequency domain. Eqs. (3.37), 
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Figure 3.2: The convergence curves of the functional J, vs. the number of propagations, for 
a state-to-state problem of a TLS. The results for the Krotov and the relaxation methods 
are presented. The number of propagations represents the amount of numerical effort. 
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Figure 3.3: The convergence curves of the functional J, vs. the number of propagations, 
for a state-to-state problem of a harmonic oscillator. The initial state is the ground state, 
and the excited state is the ground state shifted one unit to the positive direction. The 
results for the Krotov and the relaxation methods are presented. 
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(3.38) are replaced by the following equations (compare Eq. (3.19)): 

i^e^^(a;) + (1 - Ky'^'ioo) 0<K <1 



-lm{x{t)\fl\ij{t)) 



Accordingly, Eq. (3.39) is replaced by the following equations: 



(3.46) 
(3.47) 
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(3.48) 
(3.49) 



The additional update of Eq. (3.36) is added to Eq. (3.40). 

A practical remark: the first guess of the field (stage 1 in the procedure) should also be 
formulated in the frequency domain; a reasonable choice will be of the same form as /e(w). 

Now, the new method of restricting the forcing field spectrum is ready for applica- 
tion. We have chosen two examples of state-to-state problems. In both cases, the field 
is restricted to be around the frequency: uj = la.u.- We use /e(co') of the form of the hat 
function from Fig. 3.1: 



f,{uj) = Asech[20(u;-l) 



A>0 



(3.50) 



where A is a positive parameter, adjusted to our needs. 

The first example is of a TLS system, with the following unperturbed Hamiltonian: 



Hr 



1 
4 



(3.51) 



The initial state is the ground state, and the target state is the excited state. Without a 
restriction on the frequency, the preferable field for this transition is of the main resonance 
frequency, with the Bohr frequency: u = 3a.u.- However, there are other, much smaller 
resonances, in the odd fractions of the Bohr frequency. The most pronounced of them is 
at a; = la.u.- This makes this problem a test for the ability of the method to restrict the 
field spectrum and to find a satisfactory restricted field, when it is known to be possible. 

The results are shown in Fig. 3.4. The convergence curve is shown, where this time, J is 
plotted vs. the number of iterations. We see that the first guess, of the form of /e(ct;), gave 
poor results. The algorithm converges very fast to a solution, and after only one iteration, 
J attains almost its maximal value. 

The resulting e{uj) is also shown in Fig. 3.4. The spectrum has a discrete character, 
due to the use of the DCT for the finite time integral in Eq. (3.4). The resolution of e{uj) 
is determined by T, and a larger T is required for a smoother curve. The envelope shape 
of the spectrum, in the form of the hat function (Eq. (3.50)), is apparent. 

The value of Jmax is of interest — it has the physical significance of the yield in the 
process: 
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Figure 3.4: The results of the state-to-state problem of the TLS, with a restricted spectrum 
field. The convergence curve is shown in the left; the resulting spectrum is shown in the 
right. 
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Figure 3.5: The Toda potential, with a simple choice of parameters; 
V{x) = exp(— x) + X — 1 

Its value at the end of the optimization process was very close to 1: Jmax = 0.99894. 

The second problem is a more challenging one — the system is a one-dimensional 
anharmonic oscillator, with the potential: 

V{x) = exp(-x) + X - 1 (3.52) 

This potential is the so called: "Toda potential" (with a simple choice of its adjustable 
parameters; see [16]). The potential is plotted in Fig. 3.5. Its general form is close to 
that of the potential of a chemical bond; the difference between the forms is that the 
Toda potential is not dissociative, i. e. at x — > oo, also V{x) — > oo. This simplifies the 
problem, because the possibility of dissociation is prevented. This potential will be used 
also in Ch. 4 (in Subsection 3.3.3, we present a solution for the possibihty of dissociation 
in dissociative potentials, like the Morse potential). 

The initial state of the Toda potential problem is the ground-state: \(po) and the target 
state is: I^Pa). The results are shown in Fig. 3.6. e{u) is very different from the one of the 
TLS problem, but the envelope shape is still apparent. The yield was: Jmax = 0.95919. A 
better yield could be achieved if we used a smaller tolerance parameter for the convergence 
of the field. 

These examples show that the new method for quantum control with a restricted forcing 
field spectrum is applicable. In the next section, we use this method for the problem of 
harmonic generation. 

3.3 New formulation for harmonic generation 

The second part of our task is the control of the spectrum of the oscillating dipole- moment, 
which determines the spectrum of the emitted field. It is treated in a similar way to the 
first part of our task. 
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Figure 3.6: The results of the state-to-state problem of the Toda potential anharmonic 
oscillator, with a restricted spectrum field. The convergence curve is shown in the left; the 
resulting spectrum is shown in the right. 
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In Subsection 3.3.1, we introduce a new Jmaxi which represents the requirements on 
the dipole. In Subsection 3.3.2, we present the full maximization problem, and the result- 
ing Euler-Lagrange equations. Subsection 3.3.3 deals with a generalization and optional 
modifications for the new functional. 



3.3.1 New maximization functional for controlling the spectrum 
of the dipole moment 

The requirements for the response of the system to the forcing field belong to J max- Hence, 
we introduce a new Jmax that represents the special requirements of the harmonic genera- 
tion problem. 

The quantity that we want to control is the dipole-moment expectation value: 

{fi){t) = m)v^\m) (3.53) 

We want that the oscillations of this quantity, in the desired spectrum, will be of maximal 
amplitudes. The function of interest is a spectral representation of (yu)(t). We use the 
cosine transform of {fi)(t): 

= C[(/x)(t)] = (mt) cosM dt (3.54) 

Here, again, we take T to be the upper limit of integration, because the problem is undefined 

for t^T. 

{p'){oj) can attain negative values; we want to maximize the absolute value of {fi){u}) 
in the desired spectrum. It is more convenient to handle with the square of this quantity. 
Hence, we define the new Jmax in the following way: 

1 /-^ „ , .^2, 



Jmax = 2 / {^)duj ff,{uj) > (3.55) 







where f^{uj) satisfies also a normalization condition: 



n 

f^{uj) dco = 1 (3.56) 



f^{u}) has the role of a filter function. As in the case of fei^^), we can either choose a 
rectangular function, for complete filtration, or another function, for smooth filtration. 
Note that the role of the filtration in the present case, is different from that of Sec. 3.1: 
we do not require that there will be no response to the forcing field outside the desired 
spectrum, represented by f\i{oj); we simply do not encourage this response. 

The conditions for an extremal in Eqs. (2.15), (2.16), result in the following equation 
for |x(^)) (instead of (2.20) or (2.26); the derivation is given in App. A): 

= -zH(t) \x{t)) - [/,(c.)(A)H] a m)) , |x(T)) = (3.57) 
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The condition on represents the natural boundary conditions of the problem, as in 

Sec. 2.2. 

The expression: 

is a function of t. Its meaning will be more apparent if we write it in the following way: 

c-'{f,{u)cmm (3.58) 

The meaning of this expression can be described as follows: We transform {fi)(t) to the 
frequency domain; then, we apply a spectral filtration to the resulting function of u, by 
multiplying by the filter function f^{uj); the result is transformed back to the time domain. 
We can give this expression the interpretation of the filtered, time dependent dipole moment 
expectation value. 

Note the similarity of Eq. (3.57) to Eq. (2.26): Here, again, we have the form of the 
inhomogeneous Schrodinger equation; The time dependent operator w{t)0{t) is replaced 
by the time dependent operator: 

C-'{f,{oo)C[{mt)]}f^ (3.59) 

We can give a similar interpretation to the two operators: The operator w(t)0{t) is inter- 
preted as the operator 0(t), weighted by w(t) in the time domain; likewise, (3.59) may be 
interpreted as the time dependent operator {fl)(t)fl, weighted by f^{uj) in the frequency 
domain. 

Jpenal IS uot iuvolvcd in the derivation of Eq. (3.57), because it does not have an explicit 
dependence on \ip(t)). It follows, that the new Jmax may be used with various kinds of Jpenai, 
with the same Euler-Lagrange equation, (3.57). For the harmonic generation problem, we 
use the Jpenai from Eq. (3.6). 



3.3.2 The full maximization problem 



We present here, for convenience, the full maximization problem for harmonic generation. 
The maximization functional is defined as: 



J — Jmax ~\~ Jpenal ~\~ J con 
J n 







J 



penal 



J ri 



-2Re 



T 



x{t) 



tp{t) ) dt 



Uu) > 



(3.60) 
(3.61) 

(3.62) 
(3.63) 



The variables are subject to the Schrodinger equation constraints in Eqs. (2.2), (2.7). 



45 



The conditions for an extremal in Eqs. (2.14)-(2.16), together with the constraints and 
the natural boundary conditions, result in the following set of Euler-Lagrange equations: 

d\m) = _,H(t) |^(t)) , 1^(0)) = 1^0) (3.64) 



dt 

-tU{t)\x{t))-C-' Uco){fi){co) \X{T))=0 (3.65) 



dt 

H(t) = Ho-AC~'[e(c^)] 

e{co) = Uuj)C [-Im(x(t) ^(t))] (3.66) 

3.3.3 Modifications to the maximization problem 

The skeleton of the method was summarised in Subsection 3.3.2. In this subsection we 
present several modifications to the optimization problem that are occasionally necessary. 

Generalization of Jmax to an arbitrary Hermitian operator 

We start with a generalization of the problem. In the formulation that was presented, the 
dipole moment operator plays two distinct roles: 

1. The dipole moment operator serves as the operator adjacent to the control field. 

2. The dipole moment operator expectation value is the subject of the control problem. 

It is not necessary that the same operator play both roles; for example, we can use x 
polarised laser field and require the maximal response of the dipole in the y direction in 
the desired spectrum. The operator jl^ will play the first role, and jly will play the second. 

The generalization of the formulation to include such a case is trivial; the operator fl 
represents the operator adjacent to the control field, and plays the first role. The second 
role is played by an arbitrary Hermitian operator O. Eq. (3.61) is replaced by: 

Jmax = ^J^ /0(C^)(0) iu)du fo{uj) > (3.67) 

where fo{^) has the same role as ffj.{uj), and is subject to the same normalization condition 
as in (3.56). Accordingly, Eq. (3.65) is replaced by: 



^^-.H(t)|x(t))-C- 



foito){6){co 



6 m)), \xiT))=0 (3.68) 



In the present work, only the case of O = /I was implemented. 
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Including an adjustable coefficient in Jmax 

The relative weight of Jpenai in J determines tlie "cost" of large fields. The relative weight 
can be controlled, in principle, by only one adjustable coefficient — the coefficient a of 
Jpenai- In practicc, it is often recommended to include also an adjustable positive coefficient 

ill Jmax • 

Jmax = 2 A /oM(o) iLo)dco A > (3.69) 

The reason is, that increasing fe{u}) (which means decreasing the relative weight of Jpenai) 
sometimes results in a numerical instability. In such a case, it is recommended to increase 
A instead. 

It is convenient to define: 

foiu) = Xfoiu) (3.70) 

as we did in Subsection 3.1.1 in the context of the forcing field filter function. 
Eq. (3.68) is modified in the following way: 

dixit)) 



dt 



-zH(t) Ixit)) - 



fo{u;)(0)iu 



0\m) (3.71) 



Preventing the system from occupying undesirable states 

When dealing with the problem of harmonic generation for a realistic chemical-bond po- 
tential (like the Morse potential), it is important to prevent the possibility of dissociation. 
Occupation of states with eigenenergies above the depth of the well means a partial disso- 
ciation; hence, the occupation of these states should be prohibited. Moreover, the states 
with energy above the depth of the well have no physical meaning in the numerical method 
of representing the spatial variable, used in the present work (see App. B). For these rea- 
sons, it is important to include in our formulation the possibility of restricting the allowed 
state-space of the quantum system. 

To deal with the problem, we should insert two changes in the optimization problem: 

1. The forbidden components of the state vector, that may increase Jmax, are encouraged 
by the recent version of Jmax- First, we should avoid this encouragement. 

2. A penalty should be put on the forbidden states. 

The first change is easily achieved, by using only the allowed components of \ip{t)) in 
Jmax- Let us define the projection operator on the allowed subspace: 

L 

= ^ \if^) {if^l (3.72) 

n=0 

where Ei is the maximal allowed eigenenergy. The "refined" expectation value of O, 
computed only from the allowed components, is defined as: 

6\ (t) = (Pam 6 Pam) = U{t) PaOPa m) (3-73) 
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Let us define the following operator: 



P.OP. 



(3.74) 



The modified Jmnn- is: 



J n 



fo{uj)(6a) {Uj)duj 



(3.75) 



Now it is apparent, that Eq. (3.71) can be used also in this case, just by replacing O with 

The idea of penalizing the functional for forbidden states has already been suggested 
in [8]. We insert into J an additional penalty term, defined as: 



J 



forb 



-7 



i/j{t) Pf i){t) ) dt 



7 > 



(3.76) 



7 is an adjustable penalty factor; P/ is the projection operator on the subspace of forbidden 
functions, where: P f = 1 — P^. 

Our experience shows that including J forb in its present form increases greatly the 
difficulty in the optimization process. If it is practical to compute the forbidden eigenstates 
we can offer a modified version to J/orb, in which this difficulty is decreased significantly: 



J forb — 



i){t) p) i){t) ) dt 



N 



p7 



n=L+l 



7n > 



(3.77) 
(3.78) 



where is the dimension of the eigenvector. In this way, we can choose 7„ to increase 
gradually with n, and to achieve a smoother filtration of undesirable states. Note that 
when 7„ is constant in n, we return to (3.76), so (3.77) is a generalization of (3.76). 

This method of penalizing J for undesirable states should not be used with too large 7„ 
values — this might cause the maximum of J to be 0, with e{t) = 0. Large 7„ values might 
also result in a numerical instability in the propagation process. This limits the efficiency 
of the method. 

After inserting the two modifications into the functional, the Euler-Lagrange equa- 
tion (3.71) is replaced by the following one: 



dixit)) 
dt 



-zH(t) \x{t)) - 



fo{uj){6a){uj) 



p] \ \m) 



(3.79) 



Prevention of the "ringing" phenomenon 

Jmax presented in (3.60) do not include a boundary term in t = T, like the one that appears 
in (2.29). This results in the natural boundary condition: |x(^)) = 0. Using a natural 
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t 



15 



Figure 3.7: The figure demonstrates the source of the ringing phenomenon in DCT. The 
original signal, defined in the time interval: [0, T], is plotted in black; the DCT of this 
signal represents the extended periodical function (in green). This function is given by 
"folding" the signal, and treating the resulting pattern, defined in [— T, T], as a single 
period of the periodical function. Discontinuities in the first derivative (marked by red 
circles) are present in the extended function (unless the first derivative at the boundaries 
of the original signal is 0). These discontinuities contain very high frequencies, which result 
in the ringing phenomenon. 

boundary condition is, in some sense, a "waste" of an opportunity for an additional control 
requirement. In this section, we introduce a boundary term, that results in a non-natural 
boundary condition for |x(^))- 

The necessity of this boundary term is numerical, and not physical — it helps in 
preventing the "ringing" phenomenon, that may appear when using the discrete-cosine- 
transform (DCT) for the finite time integral: 



The extended periodical function, represented by the DCT of the signal, has disconti- 
nuities in the first derivative, unless the first derivative of the signal is at the boundaries 
of the signal (see Fig. 3.7). This discontinuity includes very high frequencies, which result 
in the ringing phenomenon throughout the transformed vector. The effect is much smaller 
than the one that appears when using a discrete- Fourier-transform (DFT), where there is 
a discontinuity in the periodical function itself (see [20]); however, the small effect becomes 
a trouble, when there is an interest in small amplitude zones in the spectrum. 

Usually, in our case, there is no problem with the boundary t = 0; the reason is that 
the initial state is typically chosen to be the ground state, and: 




(3.80) 








dt 
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However, there is a problem for the t = T boundary. Frequently, the harmonic generation 
effect is very small, and the noise resulting by the ringing is a serious trouble. The situation 
is even worse: The algorithm tends to increase the ringing, in order to increase Jmaxi instead 
of using a physical mechanism; the intensity of the resulting fields is greatly increased close 
to T; as a result, the dipole oscillates in a wild manner, in a way that the first derivative 
at T is maximized, and the discontinuity in the extended periodic function is increased. In 
this situation, the method completely fails. 

We propose to solve this problem by trying to force: 

d(6){T) 



dt 

This is achieved by the insertion of the following boundary term into the functional: 



J} 



bound 



d{6){T) 



dt 



K > 



(3.81) 



K is a positive parameter, that determines the relative importance of Jbound in the functional. 
This boundary term is maximized when the magnitude of the derivative at the boundary 
is minimized. Actually, it is a penalty term, and k has the role of a penalty factor. 
Taking the expectation value of both sides of the Heisenberg equation, we have: 



d(6){T) 



dt 



H(T),6 )(T) 



(3.82) 



Using Eq. (3.82), we are able to derive the following boundary condition (see App. A): 



\X{T)) 



H(T),6 )(r) H(T),6 \i,{T)) 



(3.83) 



The insertion of this Jbound into the functional might increase the difficulty in the 
optimization process. Hence, it should be used only when necessary. 



The full general maximization problem 

For convenience, we summarise all the mentioned modifications, by presenting the most 
general version of the optimization problem. 
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The maximization functional is defined as: 



J — Jmax ~l~ Jbound ~l~ Jforb ~l~ Jpenal ~l~ Jcon 

~ — ^2 





Jbound — ' 

Jforb = 
Jpenal — 
J rnn ^ 2I^e 







dt 



-e^(ci;) du 



T 



xit) 



'ip{t) ) dt 



foico) > 



K>0 



(3.84) 
(3.85) 

(3.86) 

(3.87) 
(3.88) 
(3.89) 



The resulting Euler-Lagrange equations are Eqs. (3.64), (3.66), together with the fol- 
lowing equation for |x(^)): 



-zH(t) \xit)) - 



fo{u)(6a){uj) 



dixit)) 
dt 



\x{t)) = k( h(t),6, )(t) h(t),6, |^(t)) 



o. - P] |^(t)) 



(3.90) 
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Chapter 4 

Results and discussion 



In the present chapter, the new method is apphed to several simple problems. Two classes 
of problems are discussed: 

1. Many level system problems, formulated in the basis of the Hq eigenstates; 

2. Anharmonic oscillator problems. 

The results are analysed and discussed in order to obtain more general conclusions on 
possible mechanisms of harmonic generation. 

When discussing mechanisms of harmonic generation, we should distinguish between 
two elements of the mechanism; in general, the mechanism is characterized by a steady 
state, in which the dipole spectrum contains the desired frequency. The steady state is 
characterized by a periodic pattern of changes in the system and the forcing field. The 
mechanism consists of the following elements: 

1. The mechanism of achieving the steady state; 

2. The mechanism of the process in the steady state. 

These two elements cannot always be assigned to two distinct stages in the process: using 
the new method, we frequently observe that the system passes during the process through 
many such "steady states"; any of these has its own character during a portion of the 
process, and produces the desired frequency; then, it is replaced by another, improved 
"steady state" . 

In order to be able to choose appropriate problems to test the new method, we need 
some previous insight into possible mechanisms of harmonic generation. We start from the 
second element of the mechanism, mentioned above. The characteristic frequencies of the 
system are the Bohr frequencies; the Bohr frequency of two eigenstates, \(prn) and \(pn), is 
defined as: 

(^mn = Em — En (4.1) 

It is certain that the system is able to emit radiation in the Bohr frequency of two states, 
when they are coupled by the dipole moment operator, and the higher energy level is 
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occupied. Hence, the problems are constructed in a way that the target frequency is one of 
the Bohr frequencies, Umn, and the two states are coupled by /}. The expected mechanism 
is based on the occupation of the higher energy level, Em- 

The first element of the mechanism must be based on fields with lower frequencies than 
Umn- Hence, we cannot use the main resonance frequency for the transition from \ipn) to 
\(Pm), which is equal to ujmn- The occupation transfer may be achieved by utilizing other 
resonance frequencies of the system, in two ways: 

1. Using a secondary resonance frequency for the transition from \ifn) to \(fm), i- e. using 
one of the odd fractions of Umn- 

u='^ / = 3,5,7,... (4.2) 

2. Using intermediate states for the occupation transfer; this will be demonstrated in 
Subsection 4.1.2. 

Our discussion on possible mechanisms will concentrate on the second element of the 
mechanisms. 

In Sec. 4.1, the results for many level systems are presented. In Sec. 4.2, the results of 
Sec. 4.1 are analysed. The analysis leads to general conclusions on the second element of 
harmonic generation mechanisms, mentioned above. In Sec. 4.3, the results for the more 
complicated anharmonic oscillator problems are presented. The conclusions from Sec. 4.2 
are used to analyse the results. In Sec. 4.4, the deficiencies and the problems in the new 
method are discussed. 

A general remark for all problems: the initial state is always chosen to be the ground 
state of the system. 



4.1 Many level systems 
4.1.1 Two level system (TLS) 

We start from the simplest problem — harmonic generation in a TLS. We try to see 
if the new method will follow a simple mechanism of harmonic generation; it consists of 
utilizing one of the secondary resonances of the system for transition to the excited state, 
and emission at the Bohr frequency, ui^q. 

The unperturbed Hamiltonian and the dipole moment operator are the same as in 
Eqs. (3.51), (3.45), respectively (see also Table 4.1). We want to maximize the emission at 
1^1,0 = 3a.u.. As in Subsection 3.2.3, the forcing field is restricted to be around u = la.u., 
the first odd fraction (1/3) of Ui^. /e(cu) is a hat function (see Fig. 3.1 for the general 
shape), and //^(ci;) is a Gaussian function (see Fig. 4.1). 

The details of the problem are summarised in Table 4.1. 

The optimization process converges rapidly to a solution (Fig. 4.2). The results are 
presented in Figs. 4.3-4.5. 
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Figure 4.1: f^{uj) of the TLS problem 



Ho 


[1 01 
[O 4J 




[%] 


l^o) 


[I] 


T 


100 




20 sech[20(a; - 1)^] 




exp[-10(a;-3)2] 




sech[20(w - 1)^] 




0.5 


tolerance 


10-3 



Table 4.1: The details of the TLS problem 
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1 2 3 4 5 

number of iterations 
Figure 4.2: The convergence curve of the TLS problem 
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Figure 4.3: The e(u;) curve of the TLS problem; the spectrum is characterized by a large 
distribution around uj = la.u.- The hat function envelope is apparent. 

The resulting forcing field (Fig. 4.3) is rather different from our expectations: instead of 
being consisted of a single frequency, u = la.u., ^{^) has a large distribution of frequencies 
around u = la.u.- The hat function envelope is apparent. 

The resulting {fl){u!) (Fig. 4.4) mainly consists of a large peak at uifi, as could be 
expected. There is also a small response in the neighbourhood of cui o- This is probably a 
sequence of the chosen shape of /^(w), which is not completely localized. 

Examining the time-picture of the system (Fig. 4.5) is very edifying. We observe that 
the system reaches rapidly an equal occupation of the two eigenstates. Once the equal 
occupation is achieved, the (/i)(t) oscillations reach their maximal amplitude: la.u., and 
the forcing field is "turned off". This apparently means that the system achieved its 
optimal state for the emission at Ui^. This is the steady state of the harmonic generation 
process. 
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Figure 4.4: The curve of the TLS problem; the dipole oscillates mainly at the 

Bohr frequency: Uifi = Sa.u.- A small response exists also in the neighbourhood of ouifi, in 
accordance with the Gaussian shape of f^{uj) (see Fig. 4.1) 



4.1.2 Three level system (3LS) 

The next example is of a 3 level system (3LS). The unperturbed Hamiltonian is: 



Hr 



1 
1.9 
3 



(4.3) 



The energy levels are almost equidistant, with an energy difference of around Iq.m/- The 
dipole moment operator is: 

[0 1 1" 

/i= 1 1 (4.4) 
[l 1 0_ 

This operator couples between all levels. 

In this problem, we examine the ability of the method to follow another mechanism 
for achieving the steady state, based on an intermediate state for occupation transfer. 
The field is restricted to be with frequencies around u = la.u.- This is the region of the 
resonance frequencies of the neighbouring levels; the resonance frequency of the levels Iv^o) 
and \(p2), i- e. a;2,o = 2a.u., is out of the allowed frequency region, /^(ci;) is a narrow Gaussian 
(see Fig. 4.6). We require that the dipole moment will oscillate at the Bohr frequency 002,0- 
ffiioj) has the same form as /e(a;). 

According to our requirements on the forcing field spectrum, the forcing field is at 
resonance for the transition between neighbouring levels, but is off resonant for the tran- 
sition between the non-neighbouring levels. In the expected mechanism, serves as an 
intermediate state for transferring occupation from |(^o) to \^2)- 



^The levels in this example, and in the example of Subsection 4.1.3, are not chosen to be exactly 
equidistant; this is because a system with equidistant levels possesses symmetry properties; these may 
limit the controllability of the system. The topic of controllability is discussed in [1, 2]. 
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Figure 4.5: The time picture of the system, for the TLS problem; at the top: The e(t) curve; 
in the middle: The {fl){t) curve; at the bottom: The occupation vs. time curves; once 
the two levels are equally occupied, the dipole-moment oscillations reach their maximal 
amplitude: la.u., and the forcing field is "turned off". 
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Figure 4.6: /^(w) of the 3LS problem 
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Figure 4.7: The convergence curve of the first 3LS problem 

The details of the problem are summarised in Table 4.2. The parameter a is chosen to 
be relatively large. After discussing the results of this problem, we will repeat the problem 
with another choice of a. 

The convergence of the optimization process is very fast (Fig. 4.7). The results are 
shown in Figs. 4.8-4.10. 

The resulting e{u) (Fig. 4.8) has a very apparent Gaussian envelope, in accordance with 
/e(a;). The frequencies are distributed around the region of the resonance frequencies of 
the neighbouring levels. 

The {fi){uj) spectrum (Fig. 4.9) mainly consists of the Bohr frequencies of the system. 
In accordance with our requirements, there is a large response at U2fi- Smaller response 
exists at the forcing field frequency region, with extrema at the Bohr frequencies of the 
neighbouring levels. 

The time-picture of the system is presented in Fig. 4.10. It seems from the occupation 
curve that the mechanism of harmonic generation is similar to the expected one; at the 
beginning of the process, there is a gradual increase in the occupation of \(pi), along with 
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Table 4.2: The details of the 3LS problem, with large a 
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Figure 4.8: The e{u) curve of the first 3LS problem; the frequencies are distributed around 
the Bohr frequencies of the neighbouring energy levels. The Gaussian envelope is apparent. 
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Figure 4.9: The {fi){uj) curve of the first 3LS problem; the spectrum mainly consists of 
the Bohr frequencies of the system. The largest component of the spectrum is at the Bohr 
frequency of the non- neighbouring levels: a;2,o = "^a.u.- Smaller extrema exist at the Bohr 
frequencies of the neighbouring levels, near u = la.u.- 



an increase in the occupation of \^P2)] then, the occupation of \ipi) falls off to 0, and is 
transferred to \(f2)- This implies that \(fi) serves as an intermediate state for occupation 
transfer. However, it is not absolutely certain that there is no direct transfer between \(fo) 
and |v92)- 

Despite the difference in the mechanism between the TLS and 3LS problems, similar 
features are observed in the time-pictures of the two systems. The 3LS system reaches an 
equal occupation of 0.5 in the states \ipo) and |v52), which their Bohr frequency is equal to 
the desired target frequency; then, the forcing field is "turned off" again. 

Now, we solve the same problem with a smaller a. The details of the problem are 
summarised in Table 4.3. 

The convergence curve is shown in Fig. 4.11. The resulting spectra of e{u) and 
do not seem to possess any new interesting features. However, the time-picture of the 
system is different from that of the previous problems: in this case, the forcing field is not 
completely turned off. The occupation time picture (Fig. 4.12) has a general similarity to 
that of the previous problem: the states \ipo) and \ip2) reach an equal occupation of 0.5, 
and the general picture of the system do not change much at the rest of the propagation. 
However, small occupation transfers between the levels continue, and we observe small 
"jumps" in the occupation of |v?i)- It seems that this state has a role also in the second 
element of the harmonic generation process. 

It is interesting to compare the value of Jmax of the two 3LS problems, with that 
computed using the following \ip(t)) sequence: 

m)) = exp (-zHot) [|^o) + 1^2)]} (4.5) 

In this \ip(t)) sequence, the eigenstates \ipo) and \ip2) have a constant occupation of 0.5. 
Using this I'lpit)) sequence, we have: Jmax = 24.0. In the first 3LS problem, we have a 
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Figure 4.10: The time picture of the system for the first 3LS problem; at the top: The e(t) 
curve; at the bottom: The occupation vs. time curves; the states \(fo) and \(f2) reach an 
equal occupation of 0.5, and the forcing field is "turned off". 
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Table 4.3: The details of the second 3LS problem, with a smaller a 
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Figure 4.11: The convergence curve for the second 3LS problem 
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Figure 4.12: The occupation vs. time curves of the system, for the second 3LS problem; the 
states \ipo) and \(p2) reach an equal occupation of 0.5, as in the previous problems; however, 
in this case, the occupation transfers between the levels continue, and the occupation of 
\ipi) is not completely during the rest of the propagation. 
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smaller value: Jmax = 21.5. This is not surprising, because the mechanism of harmonic 
generation seems to be based on achieving the steady state of an equally occupied quantum 
state; this is not achieved immediately at the beginning of the process. In the second 3LS 
problem, we have: Jmax = 24.3, which is even larger than that of the sequence of Eq. (4.5). 
This means that the steady state of equal occupation of 0.5 is not necessarily the one that 
yields the maximal emission, for a system with more than two states. 



4.1.3 Eleven level system (IILS) 

The next problem is a more challenging one: We want to observe the lO'th harmonic, 
using an 11 level system (IILS). This problem is based on the same principle as that of 
the 3LS problem. The purpose of this example is to test the ability of the new method in 
generating higher harmonics. 

The unperturbed Hamiltonian is: 



Hn 



2.1 



3.9 







6.1 







9.9 



11 



(4.6) 



The dipole moment operator is: 



1 

1 







1 

1 



1 

1 

1 



1 

1 

1 

1 







1 



1 



1 



1 
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(4.7) 



This fi couples between neighbouring eigenstates, and between the outer eigenstates: \ipo) 
and Iv^io)- The forcing field is restricted not to exceed the region of the resonance frequen- 
cies of the neighbouring levels: Wn+i.n = ^a.u.- We require an emission in the neighbourhood 
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tolerance 
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Table 4.4: The details of the IILS problem 

of the Bohr frequency of the outer energy levels: Wio,o = lOa.u.- We want to see if the new 
method succeeds to find the route of "climbing up" to Iv^io), in order to generate oscilla- 
tions in the requested frequency. This time, we use rectangular functions for /^(a;) and 

The details of the problem are summarized in Table 4.4. 

The convergence curve is shown in Fig. 4.13. This time, the starting point of the first 
guess is much less successful than in the previous examples. Nevertheless, the new method 
succeeds in the task of finding a satisfactory solution for the problem. 

The e{uj) curve is shown in Fig. 4.14. The main components of the spectrum are in 
the region of the Bohr frequencies of the neighbouring levels. However, there are also 
important components in other frequencies. The spectrum is cut sharply at the maximal 
allowed frequency, due to the rectangular feioj). 

The resulting {fi){uj) curve is shown in Fig. 4.15. It does not seem to contain any new 
interesting features. 

The time picture of the system is presented in Fig. 4.16. The system does not achieve 
the 0.5 occupation steady state, but gets close to it. We observed that it gets closer to this 
steady state after a larger number of iterations. It seems that the equally occupied state 
is more difficult to be achieved in this more complex system. 

4.2 Mathematical analysis of the results 

We want to get more insight into the mechanism that produces the requested frequencies 
at the steady states that were observed in Sec. 4.1. For this purpose, we analyse the 
spectrum of the expectation value of an arbitrary Hermitian operator. Then, we will be 
able to explain the observations of the previous section. The analysis will be very helpful 
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Figure 4.13: The convergence curve of the IILS problem 
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Figure 4.14: The e(co') curve of the IILS problem; the main components of the spectrum 
are in the region of the Bohr frequencies of the neighbouring energy levels, near u = la.u.- 
There are also important components in other frequencies. 
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Figure 4.15: The {fi){uj) curve of the IILS problem; the largest component of the spec- 
trum is at cJio.o = lOa.M.- Smaller components exist at the region Bohr frequencies of the 
neighbouring levels, near u = la.u.- 



for understanding the much more complex results of Sec. 4.3. 

Consider a quantum system of dimension N; first, as a matter of convenience, let us 
extend the definition of \fn), to include n<Oorn>A^ — 1: 



the ra'th eigenstate of Hq 




< n < A^- 1 
otherwise 



(4.8) 



An arbitrary operator (not necessarily Hermitian) 6 that operates in the N dimensional 
space is defined by the following set of equations: 



0<J<N-1 



(4.9) 



i=0 



(Note that using this definition, the index of the matrix elements Ojj starts from 0.) 

Now we are going to decompose 6 into a sum of operators; this decomposition will be 
helpful for the analysis of the spectrum of a Hermitian operator expectation value. 

Let us introduce a set of 2N — 1 operators: 



-(A^- 1) < n < A^- 1 



The operator q^'^ is defined by the following set of A^ equations: 



0<J<N-1 



(4.10) 



These operators will be found to be very useful for our analysis. The matrix elements of 



q*^*^ are: 



(4.11) 
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Figure 4.16: The time picture of the system, for the IILS problem; At the top: The {fi){t) 
curve; at the bottom: The occupation vs. time curves of \ipo) and |v?io); this time, the 
system only gets close to 0.5 occupation. 
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For instance, the matrix representation of q'-^-' is: 





1 
1 
1 







(4,12) 



.0 

Using Eq. (4.10), Eq. (4.9) may be written in the following way: 



N-j-1 



(4.13) 



i=0 



i=-j 



Let us define the following set of constants, characteristic to the operator 6: 



{-i<j<N-i-l)f]{0<j<N-l) 
otherwise 



(4.14) 



Let us introduce the following set of 2N — 1 operators, characteristic to the operator 6: 

d("), -{N-l)<n<N-l 
The operator do ■* is defined by the following set of N equations: 

d» 1^,) = d° 1^,) < J < iV - 1 (4.15) 

The di"^ are diagonal operators; The main diagonal of di* "* contains the i'th diagonal of 6. 
Using the operators di""*, we can write Eq. (4.13) in the following way: 



N-l 



^\v^)= E ^".q^^Mv'.) = E q^^^d^iy^.) (4.16) 

i=-(JV-l) j=_(Ar_i) 



Hence, we can write: 



N~l 



6= Yl ^^'^d 



(4.17) 



i=-(N-l) 



The q^^^di"''' operators have an important significance; The matrix elements of the operator 
q'^'^do are: 



kj 



q^'^dW ip\ = d°-SkA+j = Oi+jjSk,i+j = OkjSk,i+j (4.18) 
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For instance, the matrix representation of c^^^^d^o^ is: 









oio 

021 (4.19) 

032 

.0 

A well known example of an operator that may be written in the form of q(^^)di*^ is a, used 
in the treatment of the quantum harmonic oscillator system, a is defined by the following 
equation: 



If we define 6 as: 



(4.20) 
(4.21) 



6 = V2m(jJoX 

(m is the mass, and cuq is the characteristic frequency of the oscillator) we have: 

Okj = Vn{6kj+i + Skj-i) (4.22) 

Hence, we can write: 

a = q(-i)di-^) (4.23) 

Other well known examples are the operators 3^ and J_, used in the treatment of angular 
momentum in quantum mechanics. 

Now, consider the case when 6 = O, where O is an arbitrary Hermitian operator. It is 
easy to show that q^'^^dp is the adjoint operator of q*^*)d[*'*- 



Hence, Eq. (4.17) may be written in the following way 

N-l 

where: 



i=0 



q»d»+fq»d« 



o 



2 n = 

1 n>0 

Let us define the following set of Hermitian operators: 



f'TI 



Now, Eq. (4.25) may be written as a sum of the O'-"'^: 

N-l 

6 = ^ 6» 

i=0 
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< n < iV- 1 



(4.24) 



(4.25) 



(4.26) 



(4.27) 



(4.28) 



The O^"'' operators have an important physical significance: the operator O^*-* couples 
between the i'th nearest neighbours of the Hq eigenstates. These operators are useful when 
the energy levels are nearly equally spaced. 

The matrix representation of O'-^'', for instance, is: 



O 



10 



O 



0%, 







10 " ^21 

O32 



01^ 



LO 



(4.29) 



There are well known examples of Hermitian operators, that contain only this component 
(i = 1) in the sum of Eq. (4.28): X and P of the harmonic oscillator system, and the 
angular momentum operators, and iy, for a system that can be represented by an 
irreducible representation of the full rotation group. 
Using (4.28), we can write: 



N-l 



o 



Let us define: 

may be written as: 



Cn = {^n\tp) 
j=0 

Let us find the expression for (^O^^^^, in the terms of the c„. Using (4.18), we have: 

N-i-l 



(4.30) 
(4.31) 
(4.32) 



fc=0 j=0 



(4.33) 



From (4.27) and (4.33), we have: 



^Re 



'N-i-l 



O 



j=0 



(4.34) 



Now, we have the necessary tools for the analysis of the spectrum of ^Oy (t). The 

spectrum of (^) is the sum of the spectra of the ^O*^"-*^ (t). It will be instructive 

to express the ^O^"^^ (t) in the terms of the components of the state in the interaction 
picture. The state in the interaction picture is: 

(4.35) 



|^7(t)) =exp zHot \ij{t)) 
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We define the components of \ipj(t)): 



(4.36) 



We have: 



c„(t) = bn(t) exp{-iEnt) 



(4.37) 



First, we consider an important special case, when H(t) = Hq. In this case, the oc- 
cupation of the \ipn) is constant, and so are the bn- Using Eqs. (4.34), (4.37), we can 
write: 



6«) it) = f Re 



'N-i-l 



j=0 



k 



^ \bi+j\ \bj\ cos(wi+jjt + (/)ij) 



i=o 

= arg {b*_^jbjd^) 



(4.38) 



According to this expression, the spectrum of ^0*^*''y (t) consists of the Bohr frequencies 

of the i'th nearest neighbouring eigenstates. The amphtude of the term with the Bohr 
frequency: oOi+jj, is determined by the occupation of the eigenstates: \(pi+j), \<^j), and the 
magnitude of dfj, which represents the couphng between these eigenstates. Of course, the 

spectrum of ^O^ (t) consists of all the Bohr frequencies which have non-zero amplitudes. 

When the energy levels are roughly equally spaced, it is convenient, both conceptually and 
practically, to group together the Bohr frequencies that belong to the spectrum of any of 

the ^6(")^ (t). 

Consider the case when we want to maximize the amplitude of only one of the Bohr 
frequencies of the system: Umn] in addition, \(pm) and \ipn) are the only pair of eigenstates 
with this Bohr frequency. We ask: what is the optimal steady state for maximizing the 
amplitude, when the occupation remains constant? It is clear from Eq. (4.38), that all the 
occupation should be at the states \ fm), \^n), i-e. : 



+ \bn\ 



(4.39) 



Let us define: 



We have: 



= cos 6 
sin 9 



- - 2 



sm{29) 
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It follows that the maximal amplitude in (4.38) is achieved when 6 = n/A, and: 



(4.40) 



i. e. the optimal steady state is the equally occupied state. 

When the occupation is not constant, bn = bn{t). In this case, we should represent bn{t) 
by its spectral components. The Fourier transform of the bn{t) is: 



bJuj) 



^2n 



bnit) exp(za;t) dt 



bn(t) may be written as the inverse Fourier transform of bn{uj): 

1 



bn{t) 



bn{oo) exp{—iujt) doo 



(4.41) 



(4.42) 



Using Eqs. (4.34), (4.37), (4.42), we can write: 

bi+ji(^)bj{uj')dfj exp [i {uji+jj + u — lu') t] duo du' 



6«)W = -Re 




N-i-1 



-y. 



bi+j{uj) b.j{uj') 



arg 



6*^^.(a;)6,(u;V^ 



We see that the spectrum of ( O*^*^ ){t) consists of the frequencies: 



cos [(wi+jj + uj — uj')t + (f>ij{uj, u')] du doj' 

(4.43) 
uj' in which 



.CO 



bj{u') \dfj\ is significant. This allows the appearance of frequencies other than 

the Bohr frequencies. Usually, we observe that the main contribution to the spectrum is 
from u, u' values in which the difference: u — oo' is rather small, compared with ooi+jj- 
Hence, most frequently, the spectrum will be mainly distributed around the Bohr frequen- 
cies. If the energy levels are roughly equally spaced, the main contribution to the spectrum 
of ^O^*^^ (t) will be concentrated in a single region of the spectrum. 

We return to the analysis of the results of Sec. 4.1. The most pronounced observation of 
Sec. 4.1 was the steady state of the equal occupation of the two levels with the appropriate 
Bohr frequency. This is explained by the discussion above on the ideal steady state with 
constant occupation. In the TLS example, where there is a contribution only from one 
pair of levels, it is certain that this constantly occupied state is the ideal one in general. 

This is not necessarily true for the 3LS problem; p, of this problem (see Eq. (4.4)) may 
be decomposed into two Hermitian components: 



(4.44) 
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The above statement for the TLS problem is true if the only contribution to (fi) (u) is from 
(/i(^))(a;); however, a contribution from is also possible. It is certain that in order 

to get contribution from (/t(^))(u;), the state \ipi) must be occupied. Moreover, the Bohr 
frequencies that characterise are around u = la.u.', in order to get response in the 

neighbourhood of u = 2a.u., the occupation should not be constant (see Eq. (4.43)). These 
requirements for contribution from (/i(^))(a;) contradicts the constantly equally occupied 
picture. 

The first 3LS example, with a relatively large a, demonstrates a more energetically 
economical use of the forcing field; the altering of the occupation, required for a contribution 
from (/i(^))(a;), requires energy, and it is preferable to use the equal occupation mechanism, 
with field. In this case, there is contribution only from {jl'-'^')) (u) . In the second 3LS 
example, the smaller a allows the use of a slightly different mechanism, in which there is 
also a small contribution from (/i(^))(a;). 

In Fig. 4.17, (/i)(w) for the second 3LS is shown, along with the spectrum of the Hermi- 
tian components. The main contribution is from (/i(2))(a;), as could be expected; however, 
we observe that there is also a small contribution from (/i(^))(a;). This contribution makes 
the Jmax value larger than the maximal possible with the equal occupation mechanism. 

This simple example illustrates the effect of the generation of other frequencies than 
the Bohr frequencies characteristic to )(a;). In this case, the effect gives only a small 
contribution to Jmax- In more complex systems we may observe that the main contribution 
to Jmax involves such effects. We also frequently observe that an important contribution 
to Jmax comes from more than one Hermitian component. Typically, the most important 
contribution comes from the (//(")) (w) that their characteristic Bohr frequencies are close 
to the non-negligible part of f^ioj). 

We introduce a useful tool for finding the relative importance of the contributions of 
the Hermitian components; let us define the following set of N functionals: 

J!^L ^ I £ fMW^'i^) 0<n<N-l (4.45) 

Of course: 

N-l 

J J{n) 
'Jmax 7^ / ^ 'Jmax 

n=0 

because of the cross terms; this set of functionals has to be considered just a useful tool 
for qualitative view. 

4.3 Anharmonic oscillators 

In this section, we test the capability of the new method when dealing with the more 
complex systems of anharmonic oscillators. We deal with anharmonic potentials with the 
general form of a chemical bond potential (the Morse potential. Subsections. 4.3.1, 4.3.2), 
or a similar form (the Toda potential. Subsection. 4.3.3). 
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Figure 4.17: The composition of {fl){u), in the terms of the (/i("'))(a;), for the second 3LS 
problem; at the top: The curves of {fi){oj) and its Hermitian components; at the bottom: 
A close up view of (/i(i))(a;); the main contribution to {fl){oj) is from {fi'^'^^) (u) . However, 
there is also a small contribution from {fi^^^){u). 
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The energy levels of an anharmonic oscillator are roughly equally spaced. Hence, the 
treatment of Sec. 4.2 is appropriate here. 

We begin with a description of the structure of the dipole operator of an anharmonic 
oscillator, with the general form of a chemical bond potential. This description is necessary 
for the discussion on the results. 

First, we describe the structure of a dipole operator of the form: /i oc X. 

The dipole operator of the harmonic oscillator system has a simple structure: Only 
neighbouring states are coupled by jl, and ft = The only active Bohr frequency in 
fl is: 0Jn+i,n = (^0, where uq denotes the characteristic frequency of the oscillator. This 
structure determines the selection rules of the harmonic oscillator. 

The dipole operators of anharmonic oscillator systems have a more complex structure; 
the largest couplings are the di^, because of the similarity to the harmonic oscillator system. 
We may also have relatively large (Iq^ values, particularly for larger n values; they represent 

the deviation of the ^ipn X ipn^ from the bottom of the well, where x = 0. The c?^„ values 

for i > 1 are usually much smaller in magnitude. We also observe the following trends for 
the <„: 

1. decays rapidly with m, for m > 1. 

2. increases with n. 

The explanation for the first trend is simple: the selection rules of the harmonic os- 
cillator originates from its symmetry properties. In anharmonic oscillator systems, there 
are deviations from this symmetry. A pronounced value for with larger m, means a 
larger deviation from the harmonic oscillator selection rules, and requires greater break of 
symmetry. 

The second trend characterizes also the d^^ values of the harmonic oscillator system, 
which increase with the square-root of n. This trend is explained by the fact that higher 
energy eigenstates are characterized by a high density of (fn{x) at larger x values. Hence, 
the coupling between higher energy states represents larger amplitude phenomena, and is 
larger in magnitude. In the context of a chemical bond potential, an additional explanation 
may be given to the second trend: The deviation of the chemical bond potential from the 
harmonic oscillator symmetry is greater for larger deviations from x = 0. 

The dipole operator of the form: /i oc X, is a good approximation for a chemical bond 
system, in the region of the equilibrium state. However, this description is inadequate 
for larger deviations from equilibrium, where changes in the charge separation have to be 
taken into account. When dealing with larger deviations, the dependence of the dipole on 
X has to be described by a more complicated dipole function, fi{x). The description of the 
fl structure given above for the linear form of fi{x), might be inadequate for other forms. 
Nevertheless, the main features are the same, at least for the lower energy levels. 

When using a non-linear functional form of the deviations from the selection rules 
of the harmonic oscillator for the linear functional form, are larger in magnitude. 
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4.3.1 The HCl molecule 



The first example is the H '^^Cl molecule. We have chosen this problem, to test the method 
for a realistic choice of parameters. However, there is no intention of giving accurate 
predictions in this example. Hence, we make use of approximations that may be somewhat 
crude. 

The coordinate of the one-dimensional oscillator is the deviation of the inter-nuclei 
distance, rn-ch from the distance at the bottom of the well (r*): 

X = th-ci - r* (4.46) 

The potential of the bond was obtained by adjusting the parameters of the Morse 
potential: 

V{x) = Do [exp(-ax) - 1]^ (4.47) 

to experimental data on HCl: The atomization energy of HCl, and the frequency of vibra- 
tion, using the IR absorption frequency for the transition to the fundamental state. We 
made a few reasonable approximations. The resulting potential is brought in Table 4.5. 
The characteristic frequency of the bottom of the well is: 

a;o = 1.35 -10-1 (4.48) 

The dipole function was obtained by adjusting experimental data to a reasonable func- 
tional form. The data is the first 4 derivatives of fi{x) at equilibrium [17]: 



dx^ 



n = 1,2,3,4 

eq 



The functional form is: 



/i(x) = aix {1 - tanh [a2{x - as)""*]} (4.49) 

This functional form is intended to represent a nearly linear form at x = 0, that decays to 
at the dissociation region of the potential. This behaviour is typical to a homolitic bond 
cleavage. We made the approximation: 




eg 



The resulting system of equations was solved using the Symbolic Math Toolbox of MAT- 
LAB. The resulting function is complex; we take its real part (see Table 4.5). 

The resulting potential and dipole function are presented in Fig. 4.18. We can see that 
fi{x) indeed decays at the dissociation region of V{x). 

In our problem, we want to maximize the emission at the neighbourhood of second 
harmonic: 

W2,o = E2-Eo = 2.54 ■ 10-2 (4.50) 
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f 0.19309 Xj X 
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n(0.015-a;) 




1 


X domain 


[-0.69407, 3.51178) 


^grid 


32 


tolerance 


10-3 



Table 4.5: The details of the HCl problem 
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Figure 4.18: The approximated potential (blue) and dipole function (red) curves, for the 
HCl molecule; /i(x) decays to at the dissociation region of V{x). 
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Figure 4.19: The convergence curve of the HCl problem 

/e(cj) and /^(cj) were chosen to be rectangular functions. We restrict the field frequency 
not to exceed much uq (see Eq. (4.48)). /^(w) is chosen in a way, that from all the Bohr 
frequencies of i.e. ujn+2,n, only 072,0 is contained in the non-zero part of f^{oj). The 
other u:n+2,n are smaller, because of an anharmonic effect. We want to see if we get an 
equal 0.5 occupation of the eigenstates: [(/Jq), ^2)- 

It was necessary to restrict the allowed eigenstates in order to prevent dissociation and 
occupation of non-physical states. We use the method described in Subsection 3.3.3. 

The details of the problem are summarised in Table 4.5. 

The resulting convergence curve is presented in Fig. 4.19. 

In Fig. 4.20, the maximal |c„(t)| during the propagation is shown for all eigenstates. 
The restriction of the allowed eigenstates is shown to be successful in this case. 

In Fig. 4.21, the e(a;) and (/i)(w) curves are presented. The time-picture is presented 
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n 

Figure 4.20: The maximal |c„(t)| during the propagation for all ra; the forbidden states are 
of n > 19. The restriction of the allowed eigenstates is successful in this case. 

in Fig. 4.22. 

The spectrum of the forcing field mainly consists of the region of Uq, as could be 
expected. This is along with a large, negative, constant field {u = 0) component. This 
means, that the mean field has a negative value. This is also apparent from the time-picture. 
There is an important positive component at the first non-zero frequency: u = tt/T. This 
term completes half a period at T. Since e{u) represents the coefficients of cosine terms, 
this means that the field tends to more negative values during the propagation. This is 
also apparent from the time-picture. These observations require an explanation. 

also attains large values at the region of uq, as could be expected. There is a 

large, positive peak at = 0. This is clearly due to the positive deviations of from 

X = 0, typical to chemical bond potentials. There is an important negative component 
at a; = tt/T. It represents the gradual increase in this deviation during the process, due 
to occupation of higher energy levels; this can be seen in the time-picture, (/i) (cu) attains 
smaller, but significant, values in the neighbourhood of the second harmonic. The smaller 
values, compared with the first harmonic, are a consequence of the much smaller couplings 
of compared with 

The system seems to achieve a more-or-less steady state around t = SOOOa.u.. This 
steady state is interrupted at the end of the propagation due to boundary effects; these 
will be discussed in Sec. 4.4. 

We see that the occupation of the eigenstate at the steady state is more or less dis- 
tributed around n = 3. This indicates that a semi-classical view is not very far from reality 
in this case. This may be verified by following the x shape of {ipix, t)p during the process, 
which has a localized character. However, the shape may be very shallow, or with two 
maxima. 

Interestingly, contrary to our expectations, there is almost no occupation in |(/?o) at the 
steady state. This means, that almost all the contribution to Jmax comes from couplings 
between states with lower Bohr frequencies. We see that there are oscillating patterns in 
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Figure 4.21: The e{u) (green) and {fi){u}) (blue) curves of the HCl problem; the spectrum 
of the field mainly consists of the region of uq, besides a large, negative, constant field 
{u = 0) component. There is an important component at w = tt/T, that represents the 
tendency of the mean field to attain more negative values during the process (see Fig. 4.22). 
{fi){uj) has a large peak at a; = 0, due to the positive deviation from x = 0, typical to 
chemical bond potentials. There is an important component at u = tt/T, that represents 
the gradual increase in this deviation during the process, due to occupation of higher 
energy levels. There is a large response in the region of Uq. Smaller components exist at 
the neighbourhood of the second harmonic. The most of the important components of the 
two spectra are out of phase. 
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Figure 4.22: The time picture of the system, for the HCl problem; at the top: The e{t) 
(red) and (green) curves; at the bottom: The occupation vs. time curves for the 6 

first eigenstates; the mean forcing field becomes more negative, as higher energy levels are 
occupied. The system seems to reach a steady-state around t = 5000a.u.- This steady state 
is interrupted at the end of the propagation, because of boundary effects. At the steady 
state, e{t) and are out-of-phase. The eigenstates are more or less distributed around 

l^Ps). There is almost no occupation in \ipo). 
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the occupations of the states. These are necessary for the production of frequencies other 
than the Bohr frequencies of the involved eigenstates. This requires that e(t) will be active 
at the steady state. We see in the time-picture, that there are large oscillations of e{t) at 
the steady state, of frequencies close to uq. 

We observe that at the steady state, e{t) is out of phase with This results in 

opposite signs of e{u) and at the relevant frequencies. If we follow a semi-classical 

view, this means that the net work transferred to the system per oscillation cycle is 0. 
This is a necessary condition for a steady state, by definition. The tt phase difference is 
characteristic to a forced, undamped oscillator, when the forcing frequency is higher than 
that of the characteristic frequency of the system. 

We see that the constant equal occupation of |(/?o) and \ip2) is not favoured in this 
case. This is explained by the fact that the coupling di^Q is rather small in magnitude, 
compared with (ig^ of larger n, as was discussed in the beginning of the section. Hence, it 
is preferable to utilize the Bohr frequencies of pairs of states with larger n, while altering the 
occupation. For comparison: At the end of our optimization process, we have: Jmax = 16.9. 
If we compute Jmax using: 

m)) = ^ exp (-zHot) (l^o) + 1^2)) (4.51) 

we have: Jmax = 1-92. 

It is possible to get more insight into the mechanism of the steady state. The large, 
negative, constant component of the forcing field, has a simple explanation; it has the 
meaning of an addition of a time- independent term: kfi{x), to V{x), where k is positive. 
fi{x) is not very far from being linear, unless there are large deviations from the bottom 
of the well. The deviations in this problem are not very large. The addition of a nearly 
linear term to V{x) narrows the potential well, because of the asymmetry of V{x) around 
X = 0. The effect is of an increase in the separation between the energy levels. 

Another view of the same effect is to analyse the system in the terms of the original 
eigenstates, | The addition of a constant term kfi to Hq has a contribution to the 
diagonal of Hq, i.e. to eigenenergies En- According to the discussion at the beginning of 
the section, the values of the diagonal elements of /i — the (Iq^, increase with n, at least 
for linear fi{x). Hence, this contribution to the diagonal increases the separation between 
the energy levels. This is verified in Fig. 4.23 for /i(x) of our problem, by plotting the d'^^ 
vs. n. We see that for the occupied levels in this problem (see Fig. 4.20) is an increasing 
function of n. This means that the addition of to Hq induces a greater separation of 
the energy levels. 

The larger space between the new energy levels causes the new Bohr frequencies to be 
larger. This is an important part of the mechanism, since the original a;„+2,n for n > 
are not included in j]i{uj)- This may be verified, by solving a similar problem, without 
the inclusion of the to = region in fe{uj)- Another option is to use jl — jl^'^^ instead of jl. 
These two options were implemented in the context of a similar problem, with the Toda 
potential. In both cases, the resulting Jmax is much smaller (although it still exceeds the 
Jmax computed with (4.51)). 
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n 

Figure 4.23: The cIq^ vs. n curve of the HCl dipole operator; for the lower levels [n < 14), 
there is a gradual increase in the diagonal terms of /}. Only the lower levels are relevant 
to this problem — see Fig. 4.20. 

The field component of w = vr/T is explained by the fact, that during the process 
of achieving the steady state, eigenstates with increasing n values are occupied. Un+2,n 
decreases with n; this requires more negative fields, to increase the separation between the 
levels. 

As we saw, the components of e{uj) with u values around uq also play a role in the 
steady state. We can offer an explanation, using a classical argument: the frequency of 
the largest component of e{u) in this region is close to cui^o! this frequency is higher than 
the other Un+i,n- \'^o) is hardly occupied, so the characteristic frequencies of the oscillator 
are lower than that of the forcing field. The picture is of a forced oscillator, where the 
forcing field causes the system to oscillate with a slightly higher frequency. The overall 
system may be viewed as a new effective oscillator, with a higher fundamental frequency. 
The frequencies of the nonlinear effects of higher harmonics are also increased. 

We did not succeed to see higher harmonics than the second harmonic in this system. 
The couplings of higher harmonics are much smaller. In addition, the possibility of using 
higher levels, with larger couplings, is restricted by the requirement of the prevention of 
dissociation. There are problems in the method, which make the production of very small 
effects a difficult task. These problems will be discussed in Sec. 4.4. 

In the next two examples, we show that it is possible to see higher harmonics, when 
the system possesses appropriate properties. 

4.3.2 HCl with fictitious fi{x) 

In this problem, we show that it is possible to see higher harmonics, when using another, 
fictitious for the HCl problem. 

The experimental yLi(x) is not very far from being linear, unless there are large positive 
deviations from a; = 0. These deviations are hard to be achieved with the restriction on 
the allowed states, and only lower energy states are occupied. In order to increase the 
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Figure 4.24: The experimental HCl dipole function, and the fictitious dipole function; the 
deviation of the fictitious function from hnearity starts closer to x = 0. 

couplings of higher harmonics, we have to choose that deviates from linearity closer 
to a; = 0. We have chosen a function of the form of Eq. (4.49): 



It is designed to be similar to the experimental function at the neighbourhood of x = 0, 
but to decay rapidly closer to x = 0. The two functions are plotted in Fig. 4.24. 

The importance of this example is that molecules with rapidly decaying dipole functions 
do exist. If we get more successful results for this problem we can conclude that such 
molecules are preferable for harmonic generation. 

In this problem, we maximize the response in the region of the third harmonic. In 
Fig. 4.25, the couplings are plotted vs. n, for the experimental and fictitious jl. We 
see that at the lower states, the couplings are larger for the fictitious function. As we have 
already said, we are restricted to the lower states, in order to prevent dissociation. 

f^ioj) is again chosen in a way that only uj^fi = 3.73 ■ 10~^ is included in the maximized 
u interval, while the other cj„+3,„ are of lower frequencies. /e(w) is identical to that of the 
previous problem. 

In this problem, and in the next one, we use (/i)(T) instead of (/ia)(T) in Eq. (3.86), 
for convenience. 

The details of the problem are summarised in Table 4.6. 
The convergence curve is shown in Fig. 4.26. 

The main features of this problem are similar to that of the previous one. The resulting 
{fi){oj) curve is shown in Fig. 4.27. We can see a small response in the region of the third 
harmonic. There are large component at the first and second harmonics. 

The occupation vs. time curves for the first eigenstates, are shown in Fig. 4.28. We see 
that the system did not achieve a real steady state at t = T. However, regions of more 




(4.52) 
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Figure 4.25: The vs. n curves of the reahstic HCl dipole operator and the fictitious 
one; at the lower states, the couphngs are larger for the fictitious function. 
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Figure 4.26: The convergence curve of the HCl problem with fictitious 
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Figure 4.27: The {fj){uj) curve of the HCl problem with fictitious the response in the 

third harmonic region is marked by a red circle. There are also large components at w = 0, 
and at the first and second harmonics. 
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Table 4.6: The details of the HCl problem with a fictitious /i(x) 
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Figure 4.28: The occupation vs. t curve of the first 6 eigenstate, for the HCl problem with 
fictitious fi{x); the system has not achieved yet a real steady state. The state is far from 
being semi-classical. We can see regions with higher occupation of pairs of states, with 
n = 3 difference. A few of these regions are marked by black circles. 

ordered patterns exist. We can observe that in these regions there are pairs of states with 
n = 3 difference, with larger occupation. A few of them are marked by black circles. This 
time, the occupation is not distributed around a single n value, and the state is far from 
being semi-classical. 

4.3.3 The Toda anharmonic oscillator 

In this problem we show that when we do not have to be concerned on the possibility of 
dissociation, it is possible to see higher harmonics. 

We use the Toda potential from Eq. (3.52) (see Fig. 3.5). We still have to put a 
restriction on the allowed eigenstates, because the higher states are physically meaningless. 
Nevertheless, there is much more freedom for occupation of higher states in this problem, 
because the overall number of states is 128, compared with 32 in the HCl problems. 

When we used a linear dipole function: /i(a;) = x, we succeeded to see the 4'th harmonic. 
The results will not be shown here. 

We present the results for the Toda potential, with a non-linear dipole function: 



as for the problem of n{x) = x. The exponential term is chosen in a way that it will not 
compete with the exponential term of the potential itself. This makes this choice more 
realistic. 




(4.53) 



This function satisfies: 
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Table 4.7: The details of the Toda potential problem 

We maximize the response at the 6'th harmonic. f^{uj) is again chosen to include only: 
cug.o = 4.82a and not the other a;„+6,n- 

The details of the problem are summarised in Table 4.7. 

The convergence was found to be very slow. We observed, that this is typical for large 
fields, when using the relaxation method. We stopped the optimization process before the 
picture of the system was stabihzed. The convergence curve is shown in Fig. 4.29. 

The relevant part of the curve is shown in Fig. 4.30. We can see also higher 

harmonics. 

When we examine the harmonic generation process more carefully we find that the 
results are rather different from our expectations; the origin of the response at the desired 
frequency region is not from the 6'th harmonic, but from higher harmonics. In Fig. 4.31 
Jmax is plotted vs. n. We see that the main contribution comes from n = 9. The second 
important contribution is from n = 10. 

The explanation for these observations becomes clear when examining the occupation 
picture. The occupation vs. time picture is extremely complicated in this case, and is not 
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Figure 4.29: The convergence curve of the Toda potential problem 



1 

^0.5 

X 

Is 

~~^0.5 
-1 

4 5 6 7 8 9 

Lo (a.u.) 

Figure 4.30: The relevant part of the {fi){uj) curve of the Toda potential problem; the 
response in the maximized region is marked by a red circle. 
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Figure 4.31: Jmlx vs. n, for the Toda potential problem; the main contribution to J, 
comes from: n = 9 (the maximum) and n = 10. 
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Figure 4.32: The occupation at t = 93.4, vs. n, for the Toda potential problem; the n's of 
the 3 largest components are written beside the maxima. The n difference between these 
components is 9. There are other series of smaller components, with n differences of 9 or 
10. 

shown here; instead, we show the occupation picture at a single time point, close to the end 
of the propagation (at T = 100a,u., we observe undesirable boundary effects). In Fig. 4.32, 
the occupation is plotted vs. n. The occupied states are of large n, where the separation 
between the energy levels is considerably smaller. In this n region, the Bohr frequencies 
<^n+9,n and Un+io,n are the closest to the maximized region in the spectrum. This may be 
seen in Fig. 4.33. 

It is possible to recognize in Fig. 4.32 several series of states with relatively large 
occupation, with n differences of 9 or 10 between the states. This shows that the main 
contribution to Jmax indeed comes from the 9'th and 10 'th harmonics, as we concluded 
from the Jmax curve. 

We can conclude from the occupation picture that the state is very far from being 
semi-classical. If we follow \iIj{x, t) p, we can see that the state is split into many separated 
entities. Indeed, the semi-classical description is definitely inappropriate in this case. 

4.4 Problems in the new method 

4.4.1 The lack of treatment of the emitted field 

One obvious deficiency of our model is the lack of treatment of the emitted field. We 
assume that when {fi){t) oscillates, the field with the corresponding spectrum is emitted, 
but we ignore the effects of this emission. 

We can point at least two problems with this approach: 

1. The emission is a dissipative process. The loss of energy affects the system, and has 
to be taken into account. The lack of treatment of dissipative processes led us to non- 
physical steady-states, with constant field, in Sec. 4.1. In the realistic steady state 
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Figure 4.33: uJn+9,n and ci;„+io,n, vs. ra; in the region of the occupied states, these Bohr 
frequencies are close to the maximized region in the spectrum. 

of these problems, there must be a permanent, continuous transfer of energy from 
the forcing field into the system, to keep the optimal occupational state unchanged. 

2. The emitted radiation may interact with the system. 
4.4.2 Boundary effects 

As we have already mentioned (Subsection 3.3.3), the approximation of the cosine trans- 
form by the DCT causes boundary problems. We encounter boundary problems with both 
e{(jj) and {fi){u). 

We start from e{u). The DCT consists of cosine terms that complete an integer number 
of half-cycles in T. The definition of e{u) using a DCT forces to end the process with: 

dt 

This problem has a similar origin to that mentioned in Subsection 3.3.3: A non-zero final 
time derivative results in a discontinuity in the first derivative of the extended periodic 
function. This is impossible if we are restricted to low frequency components. 

We often observe that this condition causes undesirable effects at the end of the prop- 
agation. 

This problem is not very disturbing because the very end of the resulting process may 
be ignored. The effect of this problem on the resulting spectra, e{oj) and {jl){u), is rather 
small. 

The problem with {fi){uj) is that mentioned in Subsection 3.3.3. The solution that we 
proposed there to the problem is not ideal. We observe that the field is not likely to just 
be shifted by an appropriate phase. Instead, undesirable effects appear at the end of the 
propagation, in order to decrease d {jl){T)/dt. This problem is not very disturbing, as the 
problem mentioned for e{u). 
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A more serious problem is that the insertion of the boundary term increases the difficulty 
in the optimization process. This problem is more disturbing when dealing with small 
harmonic generation effects. This will be discussed in Subsection 4.4.3. 



4.4.3 Difficulties when the effect is small 

When the harmonic generation effect is small, the problem becomes a difficult task. When 
the ffist guess gives a negative J, and the maximal possible effect is very small, the opti- 
mization process often converges very rapidly to the solution: e{t) = 0, with: J = 0. 

In order to solve the problem, we have to try to make J positive for the ffist guess. 
The term Jmax is always positive. The penalty terms: Jpenah Jforb and Jbound are always 
negative. We can increase J in two ways: 

1. By increasing Jrr^x] 

2. By decreasing the magnitude of the penalty terms. 

We may suggest to play with the constant coefficients of the penalty terms (altering 
also the A coefficient of Jmax is not helpful for this purpose, because we are interested in 
the relative weight of the terms). Here, we have to distinguish between Jpenah and Jfarh 
or Jbound- The coefficient a of Jpenai may be decreased, if necessary. It is impossible to 
decrease 7„ or k (the coefficients of Jforb and Jbound), because it will make these penalty 
terms ineffective. Jbound is always necessary when the effect is small, as was discussed in 
Subsection 3.3.3. Jforb is also frequently necessary for small effect problems; the reason is, 
that in order to get a significant value of Jmax, we often need large fields. This may lead 
to the occupation of the forbidden states. 

In all the examples mentioned in this thesis, we used a rather arbitrary guess. The 
difficulty may be solved by using a more clever guess. It may be necessary to solve a 
former control problem, to produce a guess with the desired properties. 

We propose the following ideas for a former control problem, that may produce a larger 
Jmax value: 

1. It is possible to solve the problem of achieving an equal occupation of 0.5 of two states 
\^m), I'Pn), that their Bohr frequency cUmn is in the maximized region. We want to 
achieve this occupational state as fast as possible. We use the regular formulation 
for time-dependent control problems (Sec. 2.2), with: 



Of course, Jpenai is of the form of Eq. 3.62. 

We implemented this idea; it is effective only when the resulting effect is significant. 
The main problem with this approach is that when the new method starts with 
this guess it tends to get stuck in this mechanism (although the results continue to 
improve). As we saw, this mechanism is not necessarily the ideal one. 



O(t) = P^(,) = |0(t))(0(t)| 




(4.54) 
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2. For the anharmonic oscillator problems, we can suggest using a guess that shifts the 
occupation to another region of n values. When higher levels are occupied there are 
larger anharmonic effects. This gives a better starting point for the problem. We 
use the regular formulation of time-dependent problems, using a projection operator 
into the subspace of Mi <n< M2: 



Jpenai IS of the form of Eq. 3.62. 

This method was tried successfully for a problem with the Toda potential. 

3. Sometimes, in the solution of a problem of a low harmonic, we see also higher har- 
monics (see Fig. 4.30). Hence, we can use a problem of a lower harmonic as a former 
problem for higher harmonics. We did not try to implement this suggestion. 

The magnitudes of Jforb and Jbound may be decreased by solving a former simpler control 
problem, with the insertion of these penalty terms. If we need only the Jbound term, its 
magnitude may be decreased easily by shifting the phase of the field forward. We did not 
try to implement these suggestions. 




(4.55) 



n=Mi 
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Chapter 5 
Conclusion 



In the present work, a new theoretical method of calculation for controlling harmonic 
generation, was developed in the framework of QOCT. 

The development of the method involved coping with the more general problem of con- 
trol requirements that are formulated in the frequency domain. It has been shown that 
it is possible and, apparently preferable, to formulate these requirements in their natural 
domain — the frequency domain. The method that was developed for the harmonic gen- 
eration problem, was generalized to other control problems with frequency requirements. 

The new formulation required the use of the relaxation optimization method, which has 
not been used yet for QOCT problems. The relaxation method was found to be successful 
for the new formulation. 

The new method was applied to harmonic generation problems in simple systems. 
It was shown that the method succeeds to deal with these relatively simple problems. 
However, difficulties were encountered when the harmonic generation effect is very small. 
These difficulties originate from noise effects. Several ways to deal with the problem were 
suggested. Nevertheless, it seems that a satisfactory solution for the problem is still missing. 

The analysis of the results led to general conclusions on possible mechanisms of har- 
monic generation. Typical mechanisms in anharmonic oscillator systems were also dis- 
cussed. 

The limits of the ability of the new method are still unknown. The question on the 
ability of the method may be divided into 3 parts: 

1. The complexity of the problem: The method has not been tested yet for a "real" 
problem, with more than one degree of freedom, and realistic complexities. It is 
unknown if the method is effective for more complex cases than the simple problems 
mentioned in this thesis. 

2. Sm^all harmonic generation effects: It is unknown how small is the effect that the 
method will be able to deal with. This depends on the existence and efficiency of a 
possible solution to the noise problem. 

3. Production of high harmonics: The maximal harmonic that was attempted to be 
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achieved using the new method is the lO'th harmonic. It is not certain that the 
method will be able to produce much higher harmonics, like in the high-harmonic 
generation mechanism. 



There are many possible directions for taking the research further. Here are a few: 

• The method should be tested for real problems. 

• A satisfactory solution for the noise problem has to be found. 

• The method should be tested for the high-harmonic generation mechanism. 

• The effects of emission, missing in the present formulation, have to be taken into 
account. 

• The possibilities of improving the ability of the method, using clever guesses, should 
be investigated. 
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Appendix A 

The full derivation of the 
Euler-Lagrange equations for the new 
formulation 



The derivation will be performed for the most general case. The maximized functional is 
defined by: 



J — Jrnax ~l~ Jbound ~l~ J forb ~l~ Jvenal ~l~ Ji 
1 



max ~r 'Jbound ~r 'J forb ~r 'Jpenal ~r 'J con 
n —, r-2 



J, 



foiuj){Oa) {u))du 



Oa = PaOPa 
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Pa = ^ IV^n) (V^nl 
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e{t) cos{u!t) dt 
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T 



x(t) 



^lj{t) ) dt 



H(t) = Ho - fie{t) = Ho - A 
The constraint equations are: 

d\m) 



e{u) cos(ajt) du 



dt 

d{m\ 



dt 



-m) \m) 
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(A. 14) ensures that: 



{m\ = mm 



Assuming this, all the computations can be performed using (A. 13) only. 
The extremum conditions are: 
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Jcon is more easily handled after integrating by parts the expression: 

dip{t) \ 



xit) 



dt I 



dt 



We obtain: 
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The expression for the LHS of (A. 15), is obtained using (A.9), (A.20), (A. 12): 

^Jpenal ^ ^Jcon 



Se{Lu) 5e{uj) 5e{u) 

^Jpenal 2 _^ ^ 



^ Jcon 

6e{u) 



=2 Re 
= - 2 Im 



-e\uj 



(A.21) 
(A.22) 



T 



x{t) 



5H(t) 



i){t) ) dt 



{x(t) lAl i^it)) cos{ujt) dt 



= -2Im{C[(x(t)|Al^W)]} 
From (A. 15), (A.21), (A.22), (A. 23), we get the following expression for e(w): 

eiu) = Uu)C[-lm{xit)\fi\i;m 



(A.23) 



(A.24) 



In order to derive the LHS of (A. 16), we first write the explicit expression 

^f J max (X 

functional of li^it)): 
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The expression for the LHS of (A. 16), is obtained using (A.25), (A. 7), (A.20): 

^Jmax ^Jforb ^Jcon 
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Using (A.16), (A.26), (A.27), (A.28), (A.29), we obtain 
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Eq. (A. 17) gives the adjoint of (A. 30): 

dixit)) 



dt 



-zH(t) Ixit)) - C~' 



O 



(A.31) 



In order to derive the expression of the LHS of (A. 18), we write (A. 6) in a more useful 
form. Taking the expectation value of both sides of the Heisenberg equation, we have: 



d{6a}{T) 



H(r),6, )(T) 



Jbound becomes: 

The LHS of (A. 18) is: 

6J 



dt 



6J, 



bound 



SJr. 



6\^{T)) 5|^(T)) 6mT)) 

^Jbound 

mfT) 

^ J r,on 



Ki^iPiT) [H(T),0„J '0(T))(^A(T)| H(T),0« 
-(X(T)I 



-K 



(A.32) 
(A.33) 

(A.34) 
(A.35) 
(A.36) 



Using (A.18), (A.34), (A.35), (A.36), we obtain: 

(x(T)i =K(^iT) I [H(r),6,] I ^(T)) (^(T)| [h(t),6„^ 

H{T),6a]){T){ij{T)\ [h(T),6,] (A.37) 
Eq. (A. 19) gives the adjoint of (A.37): 

\X{T)) = K ( [h(T), 6,] ) (T) [h(T), 6,] \i,{T)) (A.38) 
Eqs. (A.30), (A.31), (A.37), (A.38), ensure that: 

{x{t)\ = [\xm^ 

Assuming this, all the computations can be performed using (A.31), (A.38) only. 

We collect the resulting equations, (A. 24), (A.31), (A.38), together with the constraint 
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(A.13): 



dt 

dixit)) 
dt 



1^(0)) = l^o) 
-zH(t) \x{t)) - I 



fo{u)(6a){u) 



PT !^ \m) 



IxiT)) = K ^ [H(T), 0„J ^ (T) [H(T), O, |^(T)) 

H(t) = Ho-/iC-i[e»] 
eH = /,MC[-Im(x(t) lAlV^W)] 

These are the Euler-Lagrange equations of the problem. 



(A.39) 

(A.40) 
(A.41) 
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Appendix B 
Numerical details 



B.l The propagator for the Schrodinger equation 

The propagator for the Schrodinger equation is based on a new, efficient and highly accurate 
algorithm, by Hillel Tal-Ezer. A new article on the propagator is to be published [18]. 



B.2 The performance of the Hamiltonian operations 

The Hamiltonian operations in the problems that depend on a spatial variable x, were 
performed using the Fourier grid method (see [19]). The operation of the x dependent terms 
in the Hamiltonian is performed in the x domain, and the operation of the p dependent 
term, i. e. the kinetic energy, is performed in the p domain. 



B.3 The choice of the x grid 

The X grids of the various problems are equidistant grids. The distance between neigh- 
bouring points in the x domain [xmin,Xmax) is: 



grid 

where Ngrid is the number of grid points. 
The X domain is chosen to satisfy: 

V{Xmin) = V{Xmax) (B.l) 

Vmax = ^ = ^ (B.2) 

2m 2m Ax^ ^ ' 

where Vmax is the maximal V{x), Pmax is the maximal momentum, and m is the mass, 

or reduced mass. (B.2) makes the maximal V{x) possible, equal to the maximal kinetic 
energy possible. 
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B.4 The approximation of the cosine transform 



The cosine transform was approximated using a discrete-cosine-transform that include the 
boundaries of the domain. It is sometimes referred as the DCT of the first kind (DCT-I, 
see [20]). The DCT-I of A''^ equidistant time points: 

ti, 1 = 0,1,. ..,Nt-l 

is defined as: 




Using this convention, DCT-I is its own inverse. In order to be consistent with the 
continuous formulation, the direct transform is multiplied by the factor: T/ y^^iVT^^Tjvr . 
The inverse transform is divided by this factor. This is also necessary for another reason: 
when using Eq. (B.3) as is, the definition of the spectral function g{u), represented by the 
transform, varies with the sampling frequency. This is corrected by the insertion of the 
factor. 

The computation of the DCT-I is performed by the FFT of the "folded" function 
(see [20]). 



B.5 The computation of J 

The gradient and relaxation methods involve the computation of the functional J, during 
the search for an appropriate parameter K. In order that these methods will be successful, 
it is important to compute the integrals in J as accurately as possible. The accuracy of 
the numerical integration should be the same as that of the propagator. 

The integration of the u dependent functions is performed like in the DCT-I. This is 
consistent with the accuracy of the u grid. 

The integration of the t dependent functions is performed by utilizing the internal 
Chebyshev time points of every time step. We describe here the integration method. 

Consider a time step with the boundaries: [t„,t„+i]. Within the time interval, there 
are Nc internal points — the boundary including Chebyshev points: 

tk = tn + ^{i-yk) (B.5) 

2/fc = cos k = 0,l,...,N,-l (B.6) 



iVc-1 



where At = tn+i — tn- 
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An arbitrary time- dependent function g{t) defined in the interval may be approximated 
by a truncated Chebyshev series: 



9{t) ~ ^rnTm{y) V = — (B.7) 

m=0 

where the Tm{y) are the Chebyshev polynomials. The am are given by: 



[N 

where hi is defined by Eq. (B.4). 

The integration of g{t) over the interval may be approximated using (B.7): 



/ dtK. — ^ ^ Orri / Tm(l/) dy = -At ^ a^Cm (B.9) 

(B.IO) 



771=0 777=0 

1 



., , m even 

777^ — 1 

m odd 



The integration of the T„i{y) is performed analytically. 
Inserting (B.8) into (B.9), we obtain: 

g{t) dt^AtJ2 Wkgitk) (B.ll) 

*" fc=0 

2 ^^-^^ 1 / krmi \ 

" (i^3T)]i E y—i) (B.12) 

The Wfc are the integration weights of the t^, for Nc Chebyshev points. They are indepen- 
dent of At. The expression for Wk is just a DCT-I of c^, within a factor, and an additional 
small difference (compare Eq. (B.3)). It may be computed efficiently, using a FFT. If all 
time steps have the same internal structure, the computation of the Wk has to be performed 
only once. 

B.6 The tolerance parameters 

The tolerance parameter of the convergence in the optimization procedures is the maximal 
allowed relative difference of the field at the end of the process. Let us denote the tolerance 
parameter of the optimization as r. For the problems with the regular Jpenai (Eq. (2.4)), 
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the convergence condition is: 



< r 



e = 



e(^o) 



For the problems with the new Jpenai (Eq- (3-6)), the condition is: 



I ^ new ^ old | 



< r 



e = 



e(wi) 



(B.13) 



(B.14) 



(B.15) 



(B.16) 



We also used a tolerance parameter for the propagator, unlike in the attached article. 
We allowed the possibility of more than one iteration in all time-steps, and not only in the 
first one (for efficiency, the average number of iterations should not exceed 1 too much; 
this may be achieved by choosing a sufficiently small At, and a sufficiently large Nc). Let 
us denote the tolerance parameter of the propagator as C,. The convergence condition is: 



U 



old I 



(B.17) 



where u is the solution vector at the edge of the time step interval. We take: C = 10 



r. 
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